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1.0  INTRODUCTION 

Motivated  by  the  increasing  demands  on  theoretical  seismology  to  under¬ 
stand  and  explain  the  propagation  of  waves  in  complex  environments, 
much  research  has  been  devoted  over  the  last  few  years  to  developing 
computationally  efficient  techniques  with  as  much  generality  as  possible. 
Although  the  differential  equations  of  motion  are  linear  in  the  field 
quantities  of  interest,  they  are  non-linear  in  terms  of  the  boundary 
conditions  for  most  realistic  structures.  This  fundamental  nonlinearity 
precludes  construction  of  the  solution  for  complex  structures  by  super- 
positon  of  the  solutions  for  simple  structures,  and  forces  one  into 
computationally  costly  schemes. 

Techniques  for  dealing  with  this  fundamental  nonlinearity  have  spanned 
the  range  from  the  crudest  classical  ray  tracing  approach  to  the 
computational -bound  finite  difference  type  methods.  However,  no  single 
technique  has  ever  proven  entirely  satisfactory  for  reasons  of  accuracy, 
completeness  of  solution,  generality  of  application,  cost  or  combinations 
thereof.  For  example,  in  cases  where  significant  diffraction  and  inter¬ 
ference  effects  require  "exact"  solutions,  finite  difference  techniques 
have  received  widespread  use.  Yet,  the  finite  difference  approach  is 
well  known  for  its  cumbersome  computational  demands  in  two  dimensions 
and  almost  insurmountable  computational  demands  in  three  dimensions 
even  on  the  fastest  computers  available  today. 

The  heavy  computation  requirements  of  the  finite  difference  type 
methods  are  created  by  the  necessity  of  refining  the  numerical  grid 
proportionately  to  the  wavelengths  of  interest  in  all  spatial  directions, 
including  regions  of  constant  material  properties.  For  problems 
involving  wave  propagation  in  irregular  layers  of  constant  material 
properties  within  each  layer,  however,  the  Boundary  Integral  Equation 
(BIE)  approach  provides  a  more  concise  and  efficient  formulation. 
Basically,  the  BIE  formulation  takes  advantage  of  the  fact  that  the 
propagation  of  waves  through  a  region  of  constant  material  properties 
can  be  treated  analytically,  leaving  only  the  interactions  at  the 
bounding  surfaces  to  be  treated  numerically.  Rather  than  imposing  a 
grid  over  the  entire  volume,  the  BIE  method  only  requires  gridding  of 
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the  interfaces  between  regions  of  constant  material  properties.  Not 
only  are  there  potential  savings  in  computational  effort  to  solve  a 
smaller  systems  of  equations,  but  the  formulation  represents  a  concise 
treatment  of  the  pertinent  physics  involved.  By  virtue  of  this 
contraction  of  information,  the  smaller  matrices  in  the  BIE  formulations 
are  much  denser  than  the  corresponding  matrices  in  the  finite 
difference  approach.  These  dense  matrices  are  typically  poorly 
conditioned  and  must  be  given  careful  consideration  during 
implementation  of  matrix  solution  techniques  to  avoid  numerical 
instabilities. 

Various  techniques  have  appeared  in  the  literature  for  dealing  with  the 
dense  matrices  in  the  BIE  approach.  One  technique  involves 
introduction  of  a  Kirchhoff  approximation  into  the  BIE  formalism  (eg., 
Berryhill,  1979;  Scott  and  Helmberger,  1982).  In  the  Kirchhoff 
approximation  the  interaction  between  neighboring  points  on  a  boundary 
is  ignored  by  locally  approximating  the  boundary  at  each  sample  point 
by  the  tangent  plane  at  that  point  and  then  using  plane  wave  reflection 
and  transmission  coefficients  to  compute  the  unknown  boundary  values. 
Sierra  Geophysics  has  added  two  important  modifications  to  this  classical 
Kirchhoff  approximation:  (1)  definition  of  a  locally  equivalent  incident 
pulse  from  the  neighboring  boundary  allowing  for  some  diffraction 
effects  which  would  otherwise  be  neglected;  (2)  formulation  in  terms  of 
a  series  of  two-layer  problems  allowing  higher  order  multiples  to  be 
computed  by  recursively  propagating  the  boundary  values  up  and  down 
the  stack  of  layers  so  that  the  computational  cost  increases 
proportionately  to  the  number  of  layers  rather  than  the  square  of  the 
number  of  layers. 

Even  with  the  improved  Kirchhoff  approximation,  one  is  still  confronted 
with  the  poor  conditioning  and  denseness  of  the  matrices  used  to 
propagate  the  boundary  values  forward  to  the  desired  positions.  Also, 
waves  involving  multiple  interactions  with  a  single  boundary  are 
neglected  except  for  the  pseudo-diffraction  effects  which  are  not  exact. 
Therefore  more  rigorous  techniques  are  required  to  deal  with  the  full 
system  of  boundary  integral  equations.  A  time  domain  treatment  of  the 
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full  system  of  equations  has  been  successfully  addressed  by  Cole  (1980) 
for  two-dimensional  acoustical  problems  in  geophysics.  Cole's  approach 
becomes  expensive  at  high  frequencies  or  for  late  arriving  signals  as 
the  product  of  the  frequency  step  times  the  time  step  must  be  less  than 
about  10  to  maintain  stability.  More  importantly,  the  formulation  does 
not  handle  the  dense  matrices  efficiently,  precluding  generalization  to 
three-dimensional  elastic  multilayered  problems.  Also,  it  is  difficult  to 
include  realistic  material  attenuation  and  to  suppress  late  arriving 
spurious  reflections  off  the  artificial  extremities  of  the  grid  using  a  time 
domain  formulation.  Significant  breakthroughs  have  been  made  at 
Sierra  Geophysics  in  deriving  an  iterative  frequency  domain  treatment 
of  the  full  system  of  equations  that  overcomes  these  shortcomings. 

The  iterative  BIE  treatment  deals  with  the  dense  singular  matrix 
equations  from  a  perturbation  point  of  view  with  respect  to  the  flat 
layer  solution.  Using  the  perturbation  approach,  all  inverse  matrix 
operations  are  reduced  to  simpler  deconvolution  operations,  which  are 
optimally  handled  using  Fourier  transform  algorithms.  Also,  spurious 
edge  effects  are  either  removed  completely  or  to  first  order  (in  which 
case  Fourier  tapering  is  used  to  remove  higher  order  effects)  depending 
on  the  slope  of  the  extremities  of  the  boundaries.  Another  feature 
that  makes  three-dimensional  problems  feasible  with  the  frequency 
domain  approach  is  the  computational  savings  experienced  in  varying 
the  numerical  grid  spacing  as  a  function  of  frequency  reducing  the  size 
of  the  matrix  equations  at  low  frequency  relative  to  high  frequency. 
Furthermore,  geometrical  spreading  and  material  attenuation  behavior 
can  be  invoked  at  high  frequency  where  appropriate  to  reduce  the  size 
of  the  matrix  equations  by  eliminating  the  noncontributing  interaction 
terms.  This  has  the  desirable  feature  of  allowing  the  iterative  BIE 
solution  to  execute  as  efficiently  as  the  approximate  Kirchhoff  or 
geometrical  ray  solutions  in  the  high  frequency  limit. 

The  basic  methodology  common  to  alt  the  BIE  techniques  is  presented  in 
Chapter  2.  The  formulation  is  completely  general  for  three-dimensional 
elastic  media  containing  any  number  of  irregularly  shaped  layers 
including  the  limit  of  zero  thickness  pinchouts.  A  discrete  grid  is 
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imposed  and  the  resulting  BIE  matrix  equations  are  discussed  in  terms 
of  denseness  and  singularity  considerations.  The  Kirchhoff  BIE 
solution  technique  developed  at  Sierra  Geophysics  as  discussed  above  is 
presented  in  Section  3.1  including  the  enhancements  to  include 
pseudo-diffraction  effects  and  permit  linear  cost  dependence  on  the 
number  of  layers.  Sample  problems  solved  with  this  Kirchhoff 
technique  are  presented  in  Section  3.2.  The  more  exact  iterative  BIE 
solution  technique  developed  at  Sierra  Geophysics  as  discussed  above  is 
presented  in  Section  4.1.  Sample  problems  solved  with  the  iterative 
BIE  technique  are  presented  in  Figure  4.2  including  comparisons 
showing  the  deficiencies  of  the  Kirchoff  technique  as  well  as  the 
convergence  behavior  of  the  iterative  approach.  A  brief  discussion  of 
the  current  status  and  future  extensions  of  the  iterative  BIE  technique 
is  discussed  in  Chapter  5  with  a  reference  list  following  in  Chapter  6. 
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2.0  GENERAL  FORMULATION  OF  BIE  METHODOLOGY 
The  boundary  integral  equations  describing  wave  propagation  through 
arbitrary  three-dimensional  elastic  multilayered  media  are  derived  in  two 
steps.  First,  the  known  characterization  of  wave  propagation  within  a 
single  irregular  layer  is  written  in  terms  of  integral  representations 
involving  the  full  space  Green's  functions  with  properties  of  that  layer. 
Second,  the  interaction  of  the  wavefield  is  simultaneously  imposed  at  all 
boundaries  to  satisfy  all  boundary  and  continuity  conditions  leading  to  a 
system  of  Fredholm  integral  equations  of  the  second  kind  for  the 
unknown  boundary  values.  Once  this  system  of  equations  is  solved  for 
the  boundary  values,  the  wavefield  may  be  propagated  from  the 
boundaries  to  all  receiver  positions  of  interest  within  a  given  layer 
using  the  integral  representations  of  the  first  step. 

The  model  geometry  for  the  wave  propagation  problem  solved  in  the  BIE 
formulation  is  depicted  in  Figure  1  by  N  irregular  layers  overlying  a 
semi-infinite  half-space.  The  layers  are  allowed  to  pinchout  but  not  to 
cross  in  this  formulation.  Each  layer  is  characterized  by  constant 
shear  and  compressional  wave  velocities  and  constant  densities. 
Material  attenuation  may  be  introduced  by  allowing  the  velocities  to  be 
complex.  Wave  propagation  within  a  given  layer  is  expressed  in  terms 
of  the  Green's  functions  for  a  full-space  with  the  properties  of  that 
layer.  The  formulation  is  not  restricted  to  constant  material  properties 
within  a  given  layer,  although  the  Green's  functions  for  that  case  are 
quite  simple.  The  formulation  is  carried  out  for  the  full  elastic  case 
and  the  corresponding  acoustic  formulation  is  obtainable  from  the 
derived  equations  by  replacing  the  vector  equations  with  scalar 
equations. 

The  first  step  in  the  formulation  is  to  write  expressions  for  the 
displacement  field  within  a  single  layer  without  consideration  of  the 
boundary  interaction.  In  layer  £  (£=1,2, . . .  ,N+1),  the  displacement 
vector  must  satisfy  the  homogeneous  (£^s)  or  inhomogeneous  (£=s) 
equations  of  motion  (depending  on  v'hether  or  not  the  source  is  located 
in  layer  £)  for  a  fu^  space  .h  properties  of  layer  £.  The 
Representation  Theorem  of  "  »stouynamics  (see,  for  example,  deHoop, 
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Figure  1 .  Cross-sectional  model  geometry  for  layered  half-space  used 
in  BIE  method  formed  by  N  irregular  layers  overlying  a 
uniform  half-space,  with  each  layer  characterized  by 
constant  material  properties.  The  source  and  receiver  can 
be  located  anywhere  in  the  medium. 
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1958)  provides  an  expression  for  the  displacement  vector  located  any¬ 
where  within  volume  V£  containing  layer  £  in  terms  of  integrals  of  the 
displacement  and  traction  field  over  the  bounding  surface  of  volume  V£ 
times  the  corresponding  Green's  functions  for  a  full-space  with 
properties  of  that  layer.  The  i-component  of  displacement  at  location 
x£  can  then  be  written  in  the  frequency  domain  using  the 
Representation  Theorem  for  a  volume  V£  bounded  by  layer  interfaces  S£ 
and  S, 


'£+r 


e*($£)  Ufo) 


/[ 


GJi(*A>i(y£)  -  hJ.(^£)u 


■(y£)]  dS(y£) 


7t 


„£  -»  x_£+l,->  x  x 

Gji(x£’y£+i)Tj  (W  •  HjicvWuj  <W 


]  dS(y£+r 


£+1 


+  6 


£s 


j  [Gji(x£’2£)fj(l£)]  dV(2£} 


(1) 


(i,j=l,2,3) 


in  which  the  summation  over  repeated  indices  is  understood,  the 
frequency  arguments  have  been  omitted  for  brevity  and 


y„ 


Gji(v*J 


H^.(x.,y  ) 

ji  £’Jm 


an  integration  point  on  bounding  surface  S^; 


the  j-component  of  the  full-space  Green’s^ 
function  displacement  vector  at  location  y 
on  surface  S  due  to  poigt  force  in  the 
i-direction  at  location  x.  with  properties  of 
layer  £; 


the  j-component  of  the  corresponding  Green's 
function  traction  vector  formed  from  the 
inner-product  of  the  kj -component  of  the  ^ 

Green's  function  stress  tensor  Gx.  at  location  y 
on  surface  S  with  the  k-compon^nt  of  the  111 

unit  upward  normal  at  point  y  (summing 
over  k=l,2,3);  m 
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Um(y  ) 
J 


£he  j -component  of  displacement  at  location 

y  on  surface  S  ; 

m’ 


i*(yj 

J  m 


=  the  j -component  of  the  corresponding  traction 
at  location  y^  on  surface  S  formed  from  the 
inner-product  of  the  kj -component  of  the 
stress  tensor  with  the  k-compon^nt  of  the 
unit  upward  normal  v"  at  point  y  (summing 
over  k=l ,2 ,3) ;  K  m 


the  j -component  of  the  source  function  at 
location  z^  anywhere  in  layer  £  (assuming  the 
source  is  a  Delta-function  in  space,  then  the 
volume  integral  reduces^to  the  evaluation  of 
the  integrand  at  point  z^); 
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£s 


0,  if  1*5 
1,  if  1— s 


,s  =  source  layer  number; 


^  ^  1,  if  inside  layer  £ 

£  (x^)  =  Jj,  if  x^  on  surface  bounding  layer  £ 

0,  if  x^outside  layer  £. 


In  Eq.  (1),  the  layer  comprising  volume  is  assumed  to  extend  to 
infinity  at  the  horizontal  extremes  to  eliminate  the  surface  integrals 
along  those  portions  of  the  surface  bounding  volume  and  the  nega¬ 
tive  sign  for  the  integral  over  surface  S«  is  associated  with  using  the 
upward  normal  v  =  -v  in  the  definition  of  the  traction  components. 

Once  the  boundary  values  for  U^y  )  and  T^y  )  are  determined  for 

j  m  j  m 

bounding  interfaces  and  S^,  Eq.  (1)  can  then  be  used  to  obtain 
the  displacement  field  at  any  point  within  layer  £.  Expressions  for 
the  full-space  Green's  functions  with  constant  material  properties  are 
given  in  Appendix  A  for  two  and  three  dimensional  wave  propagation  in 
elastic  as  well  as  acoustic  media.  This  completes  the  propagation  step 
of  the  B1E  formulation  and  what  remains  is  to  impose  the  boundary 
interaction  coupling. 


The  boundary  interaction  coupling  requires  simultaneous  satisfaction  of 
a  tractionless  free  surface  (interface  1)  and  continuous  displacements 
and  tractions  across  each  layer  interface  (2,3, . . .  ,N+1).  The  coupled 
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boundary  integral  equations  arising  from  the  zero  tractions  conditions 
along  the  free  surface  are  obtained  by  evaluating  Eq.  (1)  in  volume 
(layer  1)  at  a  discrete  number  (q-)  of  observation  points  x.  along 
surface  and  imposing  the  zero  traction  condition  T.(y^)  =  0  (j=1,2,3) 
for  all  quadrature  points  y^  on  surface  S1  to  yield: 

1  =  -  j  .?!>“}(?!)]  dS(?]) 


-J  [G]i<;r?2)T^2>  -  dS'?2> 


('.j-1,2,3). 


in  which  the  direct  source  contribution  (if  s=1)  is 

Fi(*i>  =  j  [^I’VVV] 


Equations  (2)  represent  a  simultaneous  set  of  3q^  Fredholm  integral 

equations  of  the  second  kind  for  the  same  number  of  unknown  dis- 

1 

placement  boundary  values  U. ,  j=1,2,3,  on  surface  S^,  which  are 
coupled  to  the  unknown  boundary  values  on  surface  Sg  through  the 
integral  over  surface  Sg-  Note  that  if  half-space  Green's  functions 
were  used  in  layer  1  instead  of  the  full-space  Green's  functions,  then 
the  terms  involving  the  traction  Green's  function,  H..,  would  not  be 
present  leading  to  a  simpler  form  of  the  free  surface  boundary  integral 
equations.  However,  analytic  expressions  for  the  half-space  Green's 
functions  are  only  available  for  a  special  class  of  two-dimensional 
problems  and  hence  will  not  be  introduced  for  layer  1  to  maintain 
generality  and  consistency  in  the  formulation. 
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The  coupled  boundary  integral  equations  arising  from  the  continuity 
conditions  across  each  layer  interface,  S£  (£=2,3, . . .  ,N+1),  are  obtained 
by  evaluating  Eq.  (1)  in  volumes  V£1  and  V£  (layers  £-1  and  £)  at  a 
discrete  number  (q£)  of  observations  x£  along  common  surface  S£  and 
imposing  the  natural  boundary  conditions  of  continuous  displacements 
and  tractions: 

Uj*'1^)  =  U®(x£)  and  Tj‘1(x£)  =  T*(x£) 
for  j=1,2,3  and  for  all  quadrature  points  y£  on  surface  S£  to  yield: 


2Ui (X£>  I  lGji  (x£’y£-l)Tj  (y£-l^  "  Hji  (x£’y£-l)U 

S.  , 


f1(y£)jdS(y£-] 


-J  [Gji1(x£’y£)Tj"1(y£)  -  HJ£1^£»y£)Uf“1(y£)]  dS(y£ 


+  6£-l,sFi(x£) 


IUi(x£)  =/[Gji  (x£*y£)Tj(y£)  ■  Hji(x£’y£)Uj(y£)j  dS(y£} 

S£ 

~J  j®jitx£’y£+l^Tj(y£+P  ”  Hji(x£,y£+l)Uj(y£+pj  dS(y£+l) 


+  6£>sF.(x£)  ,  (i,j=l,2,3),  (£=2,3 . N+l).  (4b) 


Equations  (4)  represent  a  simultaneous  set  of  6q£  Fredholm  integral 

equations  of  the  second  kind  for  the  same  number  of  unknown  dis- 

£  £ 

placement  and  traction  boundary  values,  U.  and  T^,  j=1,2,3,  on  surface 
S£,  which  are  coupled  to  the  unknown  boundary  values  on  surface  S£-1 
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and  S£+^  through  the  integrals  over  surfaces  _1  and  S^, 
respectively.  Note  that  when  £= 2  in  Eqs.  (4a),  then  the  term  involving 
Tj  (y^)  is  identically  zero  because  of  the  tractionless  free  surface 
conditions.  Also,  note  that  when  £=N+1  in  Eqs.  (4b),  then  the  integral 
over  surface  vanishes  by  virtue  of  the  radiation  conditions  implicit 

in  the  Green's  functions  for  the  underlying  semi-infinite  space. 
Boundary  integral  equations  (2)  and  (4)  are  evaluated  at  a  discrete  set 
of  q£  example  points  x£  on  each  boundary  S£  with  sample  points  from 
one  boundary  becoming  Green's  function  quadrature  points  for  the  next 
boundary.  When  completely  discretized,  equations  (2)  and  (4) 
represent  a  coupled  system  of  singular  Fredholm  integral  equations  of 
the  second  kind  for  the  unknown  boundary  values  along  all  layer 
boundaries.  The  singularities  occur  in  the  Green's  functions  when 
quadrature  point  ym  approaches  observation  point  xm  in  the  second 
integral  in  Eqs.  (4a)  and  the  first  integral  in  Eqs.  (2)  and  (4b). 
Singularities  can  also  occur  in  the  Green's  functions  when  two  adjacent 
layer  boundaries  intersect  and  would  be  evidenced  in  the  first  integral 
in  Eqs.  (4a)  or  the  second  integral  of  Eqs.  (2)  or  (4b),  depending  on 
where  the  intersections  occur.  An  equivalent  BIE  formulation  using 
integral  representations  for  the  traction  components  as  well  as  the 
displacement  components  leads  to  a  system  of  Fredholm  integral 
equations  of  the  first  kind  instead  of  the  second  kind.  At  first  this 
appears  to  be  a  simpler  formulation,  but  in  fact  requires  more  book¬ 
keeping  and  computational  effort  in  dealing  with  the  higher  order 
singularities  arising  from  the  derivations  of  the  stress  Green's  function 
components.  The  singularities  are  of  the  lowest  order  possible  using 
the  present  BIE  formulation  in  terms  of  Fredholm  integral  equations  of 
the  second  kind  and  are  rigorously  handled  using  the  iterative  BIE 
solution  technique  discussed  in  Section  4.1. 

To  complete  the  general  formulation,  the  discretized  boundary  integral 
equations  are  cast  into  matrix  form: 

tl] {u}  =  [H] {U}  -  [G] {T}  +  {F}  .  (5) 
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The  vectors  {U}  and  {T}  contain  all  the  unknown  boundary  values 
along  all  layer  boundaries: 

fU}  =  {U1h(U2}...,{U|H.1}  ;  {T}  =  {T2},{T3}...,{Tn+1} 

in  which  {U^}  are  subvectors  containing  the  3q^  unknown  displacement 
boundary  values  along  interface  £  (£=1,2, . . .  ,N+1)  and  {T^}  are  sub¬ 
vectors  containing  the  3q^  unknown  traction  boundary  values  along 
interface  £  (£=2,3, . . .  ,N+1).  The  vector  {F}  contains  the  direct  source 
contributions  along  interfaces  s  and  s-1  above  and  below  source  layer 
s: 

{F}  =  {0},{0},...f{Fs},{Fs+1},{0},{0},... 

in  which  {Fs}  is  the  subvector  containing  the  direct  upgoing  source 
field  evaluated  at  3qs  sample  points  in  Eqs.  (2)  or  (4b)  depending  on 
whether  s=1  or  s=2,3, . . .  ,N+1 ,  respectively,  and  {FS+1J  is  the  sub¬ 
vector  containing  the  direct  downgoing  source  field  evaluated  at  3qs+1 
sample  points  in  Eqs.  (4a)  if  s=1,2,...,N. 

The  matrix  [I]  contains  the  factors  of  (1/2)  from  the  left-hand 
sides  of  Eqs.  (2)  and  (4)  and  is  of  the  following  block  diagonal  form: 


in  which  [l^]  are  3q^x3q^  identity  submatrices  (£=1 ,2, . . .  ,N+1 );  all 
other  submatrices  of  [I]  are  zero.  The  matrix  [H]  contains  the  traction 
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Green's  function  components  times  the  quadrature  coefficients  from  the 
discretized  integrals  in  Eqs.  (2)  and  (4)  and  is  of  the  following  block 
tri -diagonal  form: 


1 

1.1 

1 

2,1 


It  »i,2> 

It  »2,21 

<3l 

[-H^H  H 


3 

3,4 


] 


[-H 


je-i 

£,£- 1 


o 


N  H£,£+l^ 


1-HN+1,n1  ^  TJ+l.N+l^ 
t_HN+l,N+l^ 


in  which  t H m )  (m=£-1  ,£,£+1)  are  3qAx3qm  submatrices  involving  the 
traction  Green's  function  components  H?.(x#  y  )  for  a  full-space  with 

J  I  A*  f  in 

properties  of  layer  n  for  n=1 ,2, . . .  ,N+1;  all  other  submatrices  of  [H] 
are  zero. 


The  matrix  (G]  contains  the  displacement  Green's  functions  times  the 
quadrature  coefficients  from  the  discretized  integrals  in  Eqs.  (2)  and 
(4)  and  is  of  the  following  block  tri-diagonal  form: 
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in  which  [G1^]  (m=£-1  ,£,£+1)  are  3q^x3qm  submatrices  involving  the 
displacement  Green's  function  components  G"j(*£»ym)  f°r  a  full-space 
with  properties  of  layer  n  for  n=1 ,2, . . .  ,N+1 ;  all  other  submatrices  of 
[G]  are  zero  and  submatrices  [G^]  (£=1,2)  are  not  used  in  the 
formulation  because  the  known  zero  tractions  along  the  free  surface 
have  been  removed  from  the  integral  equations. 

The  matrix  equation  represented  by  Eq.  (5)  must  be  solved  at  each 
frequency  value  of  interest  for  the  unknown  displacement  and  traction 
boundary  values  {U}  and  {T}  along  all  discrete  sample  points.  Once 
solved,  the  displacement  field  is  easily  obtained  at  additional  sample 
point  (if  desired)  by  evaluation  Eq.  (1)  at  locations  x^  in  layer  £  in 
terms  of  the  known  displacement  and  traction  boundary  values  from  the 
surrounding  layer  interfaces.  If  desired,  time  histories  are  obtained 
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through  discrete  Fourier  synthesis  by  storing  the  solved  boundary 
values  at  the  required  frequency  values. 

Equation  5  is  the  so-called  dense  matrix  equation  referred  to  in 

Chapter  1  and  being  of  lower  dimensionality,  but  higher  density  than 

the  corresponding  matrix  equation  obtained  from  full  volume  gridding  in 

the  finite  difference  of  finite  element  approaches.  For  purposes  of 

comparing  the  respective  storage  requirements  only,  assume  that  the 

number  of  sample  points  q^  along  all  interfaces  (£=1 ,2, . . .  ,N+1 )  is  a 

constant  equal  to  q.  Then  the  BIE  formulation  represents  a  system  of 

(6N+3)q  equations  for  the  (6N+3)q  unknown  boundary  values  in  contrast 

2 

to  being  proportional  to  q  in  finite  difference  methods  for  2-D 

calculations,  indicating  a  potentially  large  savings  in  computational 

effort  with  all  other  cost  parameters  held  constant.  This  potential 

savings  is  experienced  when  the  number  of  layers  N  is  less  than  the 

number  of  vertical  wavelengths  being  resolved  in  the  finite  difference 

approach,  which  is  almost  always  quite  substantial  especially  at  high 

frequencies  and/or  for  deep  layered  structures.  The  actual  storage 

requirements  for  the  two  techniques  are  comparable  with  the  BIE 

2 

formulation  having  (72N-18)q  +(6N+3)q  nonzero  matrix  entries  and  the 

2 

finite  difference  formulation  also  being  proportional  to  q  for  2-D 
calculations.  For  3-D  calculations,  both  methods  add  another 
multiplicative  factor  of  q  to  the  number  of  equations  and  storage 
requirements,  indicating  the  same  amount  of  potential  savings  in  3-D  as 
in  2-D  for  the  BIE  method.  Therefore  all  that  is  entailed  to  reach  the 
full  potential  in  computational  savings  is  being  able  to  efficiently  and 
stably  deal  with  the  smaller,  but  denser,  matrices  in  the  BIE 
formulation. 

As  mentioned  previously,  matrix  Eq.  (5)  has  singularities  arising  from 
the  Green's  functions  when  source  and  receiver  point  coincide  in  this 
formulation.  Therefore  careful  treatment  of  this  dense  matrix  equation 
is  required  to  avoid  potential  numerical  instabilities.  A  modified 
Kirchhoff  solution  technique  is  derived  in  Section  3.1  for  dealing  with 
this  dense  singular  matrix  equations  in  an  appropriate  manner.  Various 
local  and/or  high  frequency  assumptions  are  made  that  decouple 
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adjacent  sample  points  on  a  given  layer  boundary  resulting  in  an 
approximate  solution  applicable  where  dynamic  effects  of  interaction  of 
the  wavefield  with  a  boundary  are  not  important.  A  more  exact 
iterative  BIE  solution  technique  is  derived  in  Section  4.1  by  analytically 
subtracting  off  the  singularities  and  solving  the  system  of  equations 
from  a  perturbation  point  of  view  with  respect  to  flat  layers.  The 
iterative  BIE  solution  technique  not  only  handles  the  singularities,  but 
includes  all  boundary  interactions  and  also  effectively  deals  with  the 
denseness  of  the  matrix  equations. 
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3.0  SOLUTION  USING  MODIFIED  KIRCHHOFF  APPROXIMATION 
In  this  chapter,  a  modified  Kirchhoff  approximation  is  used  to  solve  the 
general  BIE  problem  formulated  in  Chapter  2.  The  solution  technique  is 
described  in  Section  3.1.  Included  is  a  discussion  of  two  important 
enhancements  that  make  Sierra  Geophysics'  technique  more  general  and 
efficient  than  the  conventional  Kirchhoff  treatments.  Also  the  major 
shortcomings  are  discussed  such  as  exclusion  of  head  waves  and 
diffraction  of  post-critical  waves  and  limitation  on  the  proximity  of  two 
adjacent  boundaries.  For  a  more  rigorous  treatment  of  the  general  BIE 
problem,  the  reader  is  referred  to  Chapter  4. 

Sample  applications  are  presented  in  Section  3.2  showing  the  versatility 
and  flexibility  of  the  Kirchhoff  algorithm.  In  the  first  example,  the 
Kirchhoff  algorithm  is  used  to  demonstrate  the  near-source  structural 
effects  on  body  waves  at  the  Nevada  Test  Site.  In  the  second  example, 
the  Kirchhoff  code  is  used  to  demonstrate  the  effects  of  basins  on 
earthquake  strong  ground  motions.  Comparisons  to  the  more  rigorous 
iterative  BIE  solution  are  made  in  Section  4.2. 


3.1  Methodology 

In  the  approximate  Kirchhoff  treatment  of  the  general  BIE  problem 
presented  in  Chapter  2,  the  interaction  of  the  propagating  wavefield 
with  the  boundaries  is  only  considered  locally  with  no  interaction 
between  neighboring  points  on  a  single  boundary.  At  each  boundary 
point,  the  interaction  of  the  incident  wave  with  the  boundary  is 
calculated  as  though  the  local  portion  of  the  incident  wavefield  were 
part  of  a  plane  wave  and  the  local  portion  of  the  interacting  boundary 
were  a  flat  plane  with  unit  normals  defined  by  the  tangent  to  the  local 
boundary.  Thereby  the  boundary  iteraction  reduces  to  simple 
convolutional  form  which  can  be  handled  analytically  allowing  the 
boundary  values  at  that  point  to  be  computed  from  the  local  direction 
and  amplitude  of  the  equivalent  incident  plane  wave  through  simple  use 
of  plane  wave  reflection  and  transmission  coefficients.  In  general, 
these  coefficients  change  as  a  function  of  frequency  and  position  on  the 
boundary.  The  definition  of  a  pointwise  equivalent  incident  wave  field 
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is  unique  in  Sierra  Geophysics'  Kirchhoff  solution  technique  and  is  an 
important  enhancement  for  including  some  diffraction  effects  which 
would  otherwise  be  neglected.  The  shortcomings  in  the  conventional 
Kirchhoff  approximation  arise  from  propagating  the  field  from  each 
discrete  point  on  the  adjacent  surface  to  the  local  boundary  point, 
calculating  the  boundary  values  at  that  point  using  the 
reflection/transmission  coefficients  and  then  summing  the  individual 
contributions  at  that  point  with  a  significant  number  of  those 
contributions  originating  from  post-critical  angles. 

A  second  feature  unique  in  Sierra  Geophysics'  Kirchhoff  solution 
technique  is  the  recursive  approach  to  the  propagation  of  the  boundary 
values  through  a  stack  of  layers.  The  first  step  in  the  recursive 
procedure  is  to  propagate  the  source  field  once  through  the  stack  of 
layers  applying  the  Kirchhoff  approximation  point  by  point  and  saving 
all  the  boundary  values.  Next,  the  propagation  direction  is  reversed 
with  the  upgoing  boundary  values  along  surface  1  propagated  back 
down  through  the  stack  of  layers  and/or  the  downgoing  boundary 
values  along  surface  N+1  propagated  back  up  through  the  stack  of 
layers.  This  provides  first  order  reflections  and  the  process  of 
recursively  cascading  up  and  down  through  the  stack  of  layers  is 
repeated,  updating  the  upgoing  and  downgoing  boundary  values  at  each 
interface  until  as  many  orders  of  multiple  reflections  are  included  as 
desired.  The  last  step  is  then  to  combine  the  updated  upgoing  and 
downgoing  boundary  values  at  all  boundary  points  from  which  the 
boundary  values  may  be  propagated  to  any  position(s)  of  interest  within 
any  layer  using  the  Green's  function  integral  representations  discussed 
in  step  one  of  the  general  BIE  formulation  in  Chapter  2.  This  second 
unique  feature  not  only  reduces  the  total  storage  requirements  because 
only  two  interfaces  need  be  considered  simultaneously,  but  also  allows 
the  computational  effort  to  be  linearly  proportional  to  both  the  desired 
order  number  of  the  multiple  reflections  and  the  number  of  layers  in 
the  stack  (in  contrast  to  the  number  of  layers  squared). 

This  Kirchhoff  solution  technique  may  be  viewed  as  a  ray  expansion  for 
BIE  methods  and  bears  a  strong  resemblance  to  geometrical  optics  (see, 
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for  example,  Hong  and  Helmberger,  1978).  In  fact,  geometrical  optics 
can  be  derived  as  high  frequency  approximation  to  the  Kirchhoff 
formulation,  indicating  that  with  proper  grid  resolution,  the  Kirchhoff 
approximation  will  provide  highly  accurate  results  in  domains  where 
geometrical  optics  is  valid  (eg.,  at  high  frequency  and  angles  away 
from  critical  where  there  is  little  communication  between  adjacent 
boundary  points).  In  domains  where  optics  breaks  down,  the  Kirchhoff 
approximation  still  provides  reasonable  results.  For  example,  optics 
gives  abrupt,  discontinuous  behavior  at  shadow  zone  boundaries, 
whereas  the  Kirchhoff  approximation  exhibits  continuous,  frequency 
dependent  behavior  expected  from  diffraction  theory.  The  Kirchhoff 
solution  is  most  reliable  in  situations  dominated  by  kinematic  effects, 
becoming  increasingly  inaccurate  (amplitudes  more  so  than  arrival  times) 
where  dynamic  effects  become  important.  These  dynamic  effects  are 
related  to  multiple  interactions  of  the  wavefield  with  a  boundary  such  as 
head  waves  and  diffraction  of  post-critical  waves).  In  such  cases 
where  dynamic  effects  are  important,  the  more  exact  iterative  BIE 
solution  technique  described  in  Section  4.1  is  recommended. 

The  approximate  Kirchhoff  solution  technique  can  be  derived  at  a  given 
frequency  from  the  general  BIE  formulation  in  Eq.  (5)  of  Chapter  2, 
although  the  more  conventional  derivation  is  in  terms  of  the  equivalent 
potential  field  BIE  formulation.  The  basis  of  the  Kirchhoff 
approximation  i»  to  find  a  set  of  assumptions  that  will  decouple  the 
interaction  of  the  wavefield  at  a  point  on  a  boundary  from  all  other 
points  on  that  boundary  and  cast  the  equations  into  convolutional  form. 
This  is  accomplished  by  assuming  that  the  boundary  is  sufficiently 
smooth  in  the  vicinity  of  a  sample  point  to  locally  allow  the  boundary 
geometry  to  be  represented  by  a  flat  plane  with  unit  normals  defined  by 
the  tangent  at  that  sample  point.  By  virtue  of  the  local  smoothness 
assumption,  the  boundary  values  and  Green's  functions  will  be  slowly 
varying  as  a  function  of  distance  from  this  sample  point.  In  particular, 
with  tangent  plane  representations  of  the  'ocal  boundary  geometry,  the 
traction  Green's  function  terms  which  involve  normal  derivatives  across 
the  boundary  may  be  discarded  at  all  sample  points  (i.e.,  the  term 
involving  the  matrix  [H]  in  Eq.  (5)  may  be  discarded).  It  is  further 
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assumed  that  the  tangent  plane  lies  sufficiently  close  to  the  global 

boundary  so  that  the  offsets  along  the  normal  direction  between  the 

tangent  plane  and  the  discrete  quadrature  points  may  be  disregarded  so 

that  the  integration  over  the  original  boundary  can  be  carried  out  in 

the  local  tangent  plane  coordinates.  This  reduces  all  vector  inner 

£-1  £ 

products  of  the  submatrices  [G  0],  m=£-1,£  and  [G  .],  m=£,£+1  with 

m  f  Xf  m  f  x 

subvectors  {T^},  £=2,3, . . . ,  N+1 ,  to  simple  convolutional  form. 

Therewith,  Eq.  (5)  can  be  rewritten  as 


[I] {U}  +  [G] {T}  =  {F} 


(6) 


in  which  each  row  of  matrix  [G]  corresponds  to  the  flat  plane  Green's 
function  operator  at  a  given  sample  point  obtained  by  rotating  and 
translating  the  original  Cartesian  coordinate  system  into  the  local 
tangent  plane  coordinates.  The  solution  of  this  system  of  equations  can 
therefore  be  handled  quite  effectively  by  analytically  applying  two- 
dimensional  Fourier  transforms  in  the  local  tangent  plane  coordinates, 
reducing  the  problem  to  simple  deconvolutional  division  operations. 

As  an  example,  the  form  of  the  deconvolutional  coefficients  will  be 
derived  for  interaction  at  a  boundary  in  3-D  structures.  The  Green's 
functions  in  the  local  coordinates  x,y  of  the  tangent  plane  are  obtained 
from  Appendix  A  by  letting  R  =  Ix-yl  =  (x2+y2)^.  Then  the  two- 
dimensional  Fourier  transform  of  the  Green's  function  G(x,y)  takes  the 
form 


g(k1(k2 


-ik^R  ikjX  ik2y 
(1/R)e  e 


dx  dy, 


(7a) 


in  which  ka  is  the  wavenumber  defined  by  frequency  (u>)  divided  by 
acoustic  velocity  (a)  of  the  layer.  Then  introducing  the  Sommerfeld 
representation 
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ikaR  C 

~ir~  J 


k  Jn(kR)dK 

L  ^ 


,.2  i2  ^ 

o  (k  -ka  V 


(7b) 


into  Eq.  (7a)  and  using  the  identity 


/' 


2  2k  ik^x 
JQ[k(x  +y  j  ]  e  dx  = 


2  cos[y(k2-k2)'5] 

\  (k2-^ 


to  integrate  over  x  leads  to: 


,  k>k 


,  0<k<k 


(7c) 


g(klfk2) 


00  < 

■// 

-®  k. 


2  cos[k2-k2)^]  k 


(k2-k^ 


ik?y 

e  dk  dy. 


(7d) 


2  2  ^ 

Finally,  introducing  the  change  of  variables  k^Ck^-k^ ) z,  the  integral 
over  y  in  Eq.  (7d)  requires  that  k^k^  for  a  nontrivial  integral  leading 
to  the  final  form  for  the  deconvolutional  coefficients  in  the  transform 
domain: 


g(krk2)  =  2n(k2+k2-k2f!*  =  -2ni/va 


(7e) 


in  which  v  is  the  vertical  wavenumber  for  acoustic  velocity  a  (i.e., 
cosine  of  the  angle  of  emergence  divided  by  a).  If  one  were  dealing 
with  potentials  instead  of  displacements  and  tractions,  this  would  allow 
the  unknown  potentials  to  be  solved  in  terms  of  the  familiar  plane  wave 
reflection  and  transmission  coefficients  at  the  artificial  tangent  plane 
boundary.  For  example,  consider  the  solution  of  the  downgoing 
potential  field  and  its  normal  derivative  d$^/dn  in  the  transform 
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domain  (k^k^w)  at  each  point  along  a  boundary  S£  with  a  source  in 
layer  £-1.  The  pair  of  boundary  integral  equations  for  the  two 
unknowns  at  this  point  takes  the  form 


WvM  * 


dWkl'k2> 

dn 


=  2F, 


.  2  Ckj ,k2) 


<C£(krk2) 


dn 


=  0 


(7  f ) 


Then  applying  the  natural  boundary  conditions 
,  d<t,£-l  P£-l  d<^£ 

V 1  "  *£  and  SE~  =  diT 


(  7  g ) 


equation  (7f)  reduces  to  two  equations  in  the  two  unknowns  $£  and 
d<j>£/dn,  which  immediately  yields  to  the  solution: 


4>£^kl’k2^  =  T£_1  £  F£-1  =  (1+R2_1  £)  F£_1 

d<}>fl 

dn"  ^kl,k2^  ~  lVa£T£-l,£  F£-l 


(7h) 


in  which 


V/p*-> '  Vp* 


2  v 


2-1, £ 


va  /P£-1 
a£-l 


+  v  /p  ’  T£-l,£ 

a£  * 


v/^'1 


+  VP* 


(7i) 


are  the  plane  wave  reflection  and  transmission  coefficients  at  boundary 
S£  given  in  terms  of  ratios  of  vertical  wavenumber  divided  by  density, 
vc/p£'  ^or  ,aYers  ^  and  respectively;  the  vertical  wavenumbers, 
va,  are  in  turn  given  in  terms  of  the  wavenumber  transform  variables 
(k^k^),  the  frequency  u)  and  the  acoustic  velocities  a£_^  and  a£, 
respectively,  as  defined  in  Eq.  (7f). 
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Returning  to  the  system  of  Kirchhoff  integral  equations  in  Eq.  (6)  for 
the  unknown  displacement  and  traction  boundary  values,  the  solution 
could  be  calculated  directly  by  solving  the  system  of  equations  in  the 
wavenumber  transform  domain  to  yield 

{U,T}  =  [I,Go)'1{F}  (8) 

in  which  the  identity  and  convolutional  submatrices  have  been  combined 
into  a  square  matrix  and  the  unknown  displacement  and  traction 
boundary  values  have  been  combined  into  a  single  vector.  The  inverse 
matrix  operation  appearing  in  Eq.  (8)  actually  reduces  to  a 
deconvolutional  problem  in  the  wavenumber  transform  domain;  however, 
the  total  computational  effort  is  proportional  to  the  number  of  layers 
squared.  To  reduce  this  squared  dependence  on  the  number  of  layers, 
the  Kirchhoff  solution  is  modified  so  that  only  two  layers  need  to  be 
considered  at  a  time  leading  to  a  linear  dependence  on  the  number  of 
layers. 

The  reduction  to  linear  cost  dependence  on  the  number  of  layers  is 
made  possible  by  realizing  the  Eq.  (6)  is  only  valid  locally  or  at  high 
frequencies  and  therefore  the  coupling  from  interfaces  separated  by 
more  than  one  layer  can  be  ignored  without  affecting  the  Kirchhoff 
solution  substantially.  To  obtain  the  first  order  Kirchhoff  solution,  the 
direct  source  field  is  propagated  once  through  the  stack  of  layers  with 
the  boundary  interaction  considered  one  interface  at  a  time  involving 
only  the  propagated  field  from  the  previous  boundary  (or  the  direct 
source  field  for  the  interface  above  or  below  the  source)  and  the 
boundary  values  along  that  interface.  If  the  source  is  located  in  the 
underlying  half-space  (layer  N+1),  this  gives  the  upgoing  boundary 
values  for  the  direct  arrivals  only;  conversely  if  the  source  is  located 
in  layer  1,  this  gives  the  direct  downgoing  boundary  values;  for 
sources  in  intermediate  layer  £,  this  gives  the  direct  upgoing  boundary 
values  along  surfaces  £,£-1,...,1  and  the  direct  downgoing  boundary 
values  along  surfaces  £+1  ,£+2, . . .  ,N+1 . 
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Assuming  that  the  source  is  located  in  layer  £,  the  direct  upgoing 
boundary  values  at  interface  Z  above  the  source  are  found  by  solving 
the  following  pair  of  equations 


2^U^0  +  *G2,£^T£^0  =  ^ 

+  <9 

in  the  wavenumber  transform  domain  at  each  decoupled  sample  point 
analogously  to  the  example  in  Eqs.  (7a)  through  (7i);  if  the  source  is 
in  layer  1,  only  the  second  equation  is  present.  The  analogous  pair  of 
equations  for  the  direct  downgoing  boundary  values  at  interface  £+1 
(£=1 ,2, . . . ,N)  are  given  by 


2*U£+1^0  +  tG£+l,£+l^T£+l^0  ^F£+i* 


(9b) 


lfjt  r/>£+l 

2‘Vl1.  +  tl£+l,£+l 1  ll£+l 


HT0A1}a  =  {0} 
0 


The  vector  subscripts  0  in  Eqs.  (9a)  and  (9b)  represent  the  zero  order 
Kirchhoff  solution  and  the  vector  superscripts  u  and  d  represent  the 
upgoing  and  downgoing  boundary  values,  respectively.  The  zero  order 
(n=0)  upgoing  boundary  values  for  Eq.  (9a)  are  propagated  up  to  the 
free  surface  by  solving  the  following  pair  of  equations  in  the  wave- 
number  transform  domain  one  interface  at  a  time  (i=£-1  ,£-2, . . .  ,2,1): 


i{U.}u  +  (G*  .HT.}U 
2l  l'n  1,1  i’n 

i{U.}U  +  [G*~*] {T. }U 
2l  i*n  i,i  i  n 


=  {0} 


(10a) 


For  i=1 ,  only  the  second  equation  in  Eq.  (10a)  is  present  and  again  the 
displacements  equal  twice  the  upgoing  field.  In  all  cases,  the  right- 
hznd  side  of  Eq.  (10a)  is  known  from  the  solution  at  the  previous 
interface.  The  analogous  pair  of  equations  for  the  zero  order  (n=0) 
downgoing  boundary  values  for  successive  interfaces  i=£+2,  £+3,...,N+1 
(£<N): 
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i{U.}d  +  [Gi_!]{T.}d 
2  i’n  1  1,1  in 

|{U.}d  +  [G1  . ] {T . }d 

2  in  l.i  in 


=  [Hi_1  .]{U.  Jd  -  [Gi_1  ]  {T.  _}d 

i,i-l  l-l  n  1,1-1'*  i-1'n 

=  {0} 


(10b) 


At  this  stage,  the  upgoing  and  downgoing  boundary  values  have  been 
calculated  and  stored  at  all  sample  points  for  the  zero  order  (n=0) 
solution,  which  represents  the  direct  arrivals  from  the  source  (no 
higher  order  multiple  reflections).  If  higher  order  multiples  are 
desired,  the  recursive  cascading  process  is  performed.  The  next  order 
solution  (n)  is  obtained  by  reversing  the  propagation  direction  from  the 
previous  pass  through  the  layers.  The  downgoing  boundary  values 
from  interface  N+1  in  Eq.  (10b)  are  used  in  the  right-hand  side  of 
Eq.  (10a)  for  i=N  and  the  upgoing  boundary  values  are  calculated 
interface  by  interface  (i=N,N-1 , . .  .2,1 )  for  the  next  order  reflection 
and  added  to  any  previously  calculated  upgoing  boundary  values. 
Analogously,  the  upgoing  boundary  values  from  interface  1  in  Eq.  (10a) 
are  used  in  the  right-hand  side  of  Eq.  (10b)  for  i=2  and  the  downgoing 
boundary  values  are  calculated  interface  by  interface  (i=2,3, . . .  ,N,N+1) 
for  this  next  order  reflection  and  added  to  any  previously  calculated 
downgoing  boundary  values. 


Note  for  the  first  order  (n=1)  reflections  only,  the  direct  arrival  terms 
(n=0)  are  zero  for  the  upgoing  boundary  values  along  surfaces  below 
the  source  (£+1 ,£+2, . . .  ,N)  and  for  the  downgoing  boundary  values 
along  surfaces  above  the  source  (£,£-1 , . . . ,1). 


After  the  cascading  process  is  complete  up  to  the  desired  reflection 
order,  the  total  boundary  values  for  the  modified  Kirchhoff  solution  are 
obtained  by  adding  the  upgoing  and  downgoing  boundary  values  along 
each  interface,  i=1 ,2, . . .  ,N+1 : 


{U.}  =  {U. }u  +  {U.}d 

*  i  *  iJn  *  l'n 

{T.}  =  {T. }U  +  {T. }d 

l  *  I’n  in 


(11) 


Then  the  displacement  field  at  any  additional  locations  within  any  layer 


SGI  -R-83-083 


28 


is  obtained  at  this  frequency  by  introducing  the  boundary  values  from 
Eq.  (11)  into  the  representation  theorem  of  Eq.  (1).  Time  domain 
results,  if  desired,  would  be  obtained  through  Fourier  synthesis  by 
calculating  modified  the  Kirchhoff  solution  at  a  discrete  set  of  frequency 
values. 

Alternatively,  one  can  formulate  the  Kirchhoff  solution  directly  in  the 
time  domain  using  time  domain  analogs  of  the  full-space  frequency 
domain  Green's  functions  given  in  Appendix  A  to  perform  the 
propagation  between  boundaries.  The  boundary  interaction  at  a  local 
point  on  surface  is  accomplished  by  performing  a  plane  wave 
decomposition  of  the  incoming  wave  field  in  the  local  tangent  plane 
coordinates.  The  effective  angle  of  incidence  of  the  equivalent  incoming 
plane  wave  at  a  local  boundary  point  is  determined  from  the  following 
plane  wave  identity  in  the  local  tangent  plane  coordinates 

d<|>£/dn  =  -  iwv^d^/dt 

by  solving  for  the  vertical  wavenumber  (i.e.,  cosine  of  angle  of 
emergence  divided  by  velocity)  at  each  time  point.  This  permits 
pointwise  application  of  the  reflection  and  transmission  coefficients  to 
the  equivalent  incoming  wave  field  from  an  adjacent  boundary,  reducing 
spurious  arrivals  from  integration  points  beyond  critical  angle. 


3.2  Results 

A  time  domain  Kirchhoff  algorithm  has  been  developed  using  the 
methodology  outlined  in  Section  3.1  and  has  been  tested  and  applied  to 
a  number  of  interesting  geologic  models.  The  algorithm  is  quite 
general,  allowing  for  three-dimensional  multilayered  elastic  model 
specifications.  The  limitations  on  the  model  input  include  a  requirement 
that  the  minimum  distance  between  two  boundaries  be  greater  than  the 
surface  sampling  distance  and  a  precaution  that  the  results  become 
unreliable  for  steeply  dipping  structures  or  any  other  situation 
involving  multiple  interactions  of  the  wavefield  with  the  boundary.  A 
direct  comparison  is  made  in  Section  4.2  between  the  Kirchhoff  solution 
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and  the  more  rigorous  iterative  BIE  solution,  graphically  displaying 
some  of  these  deficiencies.  The  Kirchhoff  algorithm  performs  all  its 
operations  for  fully  3-D  models.  However,  the  demonstrational 
applications  presented  in  this  section  involve  models  with  constant 
material  properties  perpendicular  to  the  vertical  2-D  plane  defining  the 
layer  properties. 

In  the  first  application,  the  Kirchhoff  algorithm  was  used  to  show  the 
near-source  structural  effects  on  body  waves  in  the  vicinity  of  Yucca 
Flats  at  the  Nevada  Test  Site  (Mellman,  et  al.,  1982).  The  simulated 
results  were  shown  to  reproduce  the  gross  features  observed  in  wave¬ 
form  variation  and  amplitude  variation  across  the  northern  portion  of 
the  basin,  solely  by  modeling  the  interference  effects  from  the  layer 
discontinuity  along  the  Cenozoic-Paleozoic  basement  contact.  In  the 
southern  portion  of  the  basin,  however,  the  steeply  dipping  sidewalls 
cast  some  doubts  on  the  reliability  of  the  Kirchhoff  results  and  the 
comparisons  between  synthetics  and  observations  were  less  satisfactory. 

A  contour  map  of  the  depth  to  the  Paleozoic  contact  at  Yucca  Flats  has 
been  determined  by  Herrin  and  Goforth  (1981)  as  shown  in  Figure  2. 
The  line  A-A1  is  chosen  to  be  perpendicular  to  the  major  axis  of  the 
basin,  in  a  region  which  is  essentially  two-dimensional  in  nature.  This 
allows  for  results  obtained  on  line  A-A'  to  be  applied  to  nearby  source 
locations  with  some  degree  of  confidence.  A  cross-section  view,  along 
line  A-A',  of  the  model  is  shown  in  Figure  3.  Squares  represent  grid 
points,  while  the  normals  to  the  surface  at  each  grid  point  designated 
by  a  short  line  segment  at  each  grid  point.  Note  that  irregular  grid 
spacing  has  been  deployed  to  provide  improved  accuracy  and  efficiency. 
Source  locations  are  chosen  along  line  A-A'  and  are  uniformly  spaced  at 
a  depth  of  550  meters  as  shown  on  the  artificial  intermediate  layer 
boundary  with  source  locations  numbered  from  west  to  east.  The  sides 
of  the  basin  are  tapered  at  a  depth  of  200  meters  to  prevent  numerical 
instabilities  associated  with  boundaries  being  too  close.  The  material 
properties  in  the  tuff  layer  represent  an  average  of  saturated  and 
unsaturated  values,  because  the  water  table  interface  occurring  at  a 
constant  depth  of  about  500  meters  was  removed  from  the  model  to 
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Figure  2.  Contour  map  of  Cenezoic-Paleozoic  contact  for  Yucca  Flats 
are  determined  by  Herrin  et  ai.  Line  A-A'  corresponds  to 
cross  section  shown  in  modeling  plots.  Locations  1-9  are 
selected  events  used  in  waveform  comparisons.  • 
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accommodate  the  basin  tapering.  The  P  and  S  wave  velocities  and 

density  in  the  tuff  layer  are  accordingly  set  to  2.4  km/sec,  1.4  km/sec 
3 

•  and  2.2  gm/cm  ,  respectively,  and  the  corresponding  values  in  the 
Paleozoic  zone  are  4.8  km/sec,  2.8  km/sec  and  2.6  gm/cm3. 

Reciprocity  arguments  are  involved  such  that  incident  planes  waves  with 
ray  parameter  of  .07  sec/km  provide  synthetic  seismograms  at  tele- 
seismic  receiver  distances  of  about  60  degress  due  to  the  source 
locations  displayed  along  line  A- A'  of  Figure  3.  The  synthetic  seis¬ 
mograms  are  shown  in  Figures  4,  5  and  6  for  receivers  located  in 
azimuths  west,  north  and  east  of  the  basin,  respectively.  In  each 
case,  the  response  at  the  far  field  station  is  displayed  due  to  every 
fifth  source  location  on  line  A- A*.  The  time  series  correspond  to 
realistic  explosion  source  seismograms  because  the  impulse  response 
Kirchhoff  solution  has  been  convolved  with  a  von  Seggern  and 
Blandford  (1972)  time  function  characterized  by  k=10  and  B=2,  an 
attenuation  operator  (Futterman,  1966)  characterized  by  t*=.6,  a  typical 
receiver  function  (Lundquist  and  Kam,  1982)  and  a  KS  instrument 
response.  In  all  cases  reflections  of  up  to  order  two  have  been 
included  by  making  three  passes  through  the  layers. 

The  large,  early  arriving  reflection  from  source  1  interferes 
constructively  with  the  multiple  reflections  from  the  western  portion  of 
the  basin  in  Figure  4,  giving  rise  to  the  focussed  (narrow)  first  trough 
and  second  peak  in  the  waveform.  The  greater  delay  of  this  reflection 
off  the  bottom  and  side  of  the  basin  for  sources  located  closer  to  the 
center  of  the  basin  (eg.,  sources  11  and  16)  results  in  defocussing 
(broadening)  of  the  waveform,  an  attendant  decrease  in  amplitude,  and 
the  development  of  an  inflection  in  the  second  peak.  Near  the  eastern 
portion  of  the  basin  (eg.,  sources  21  and  26),  there  is  an  increase  in 
amplitude  and  decrease  in  pulse  width  due  to  the  proximity  of  the 
eastern  wall.  The  corresponding  synthetic  seismograms  for  the  north¬ 
ern  and  eastern  azimuths  are  shown  in  Figures  5  and  6.  In  general 
there  is  a  somewhat  reduced  amplitude  and  waveform  variation  in  these 
azimuths  as  a  function  of  source  location  across  the  basin.  The 
predicted  amplitude  and  waveform  variations  in  all  three  azimuths  are 
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NTS  BASIN-WEST  AZIMUTH 


EXPLOSION  SOURCE 


Figure  4.  Synthetic  seismograms  in  the  Western  azimuth  for  model  in 
Figure  3  using  SRO  instrument,  von  Seggern-Blandford 
source  function  with  k=10,  B=2,  receiver  function  for 

station  E  from  E.  Kazakh  events  and  additional  attenuation 
operator  with  t*=.6  seconds. 
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consistent  with  observations  from  the  nine  underground  nuclear 
explosions  located  near  line  A-A1  (refer  to  Figure  2  for  the  locations 
and  to  Mellman,  et  al.,  1982,  for  the  comparisons  between  observations 
and  synthetics).  The  synthetic  seismograms  located  in  a  southern 
azimuth  (not  shown)  were  found  to  be  less  satisfactory  in  reproducing 
the  genera!  trends  in  the  observations  in  that  region.  This  is 
presumably  due  to  the  presence  of  extreme  dips  in  the  near-source 
structure  of  that  area,  which  is  not  properly  handled  with  the 
Kirchhoff  solution  technique.  The  Kirchhoff  code  proved  instrumental 
in  relating  the  observed  amplitude  and  waveform  variations  at  Yucca 
Flats  to  interference  effects  of  known  shallow  structure  rather  than  to 
focussing  caused  by  deeper  structure. 

In  a  similar  study,  the  effects  of  basins  on  earthquake  strong  ground 
motions  were  investigated  using  the  Kirchhoff  algorithm  (Hadley,  1982). 
Figure  7  shows  an  example  calculation  for  a  simple  basin  model  with  the 
earthquake  source  located  beneath  the  basin.  The  incident  energy  is 
from  the  right  side  of  the  basin,  propagating  parallel  to  the  indicated 
arrow.  The  P-wave  velocity  in  the  layer  is  2  km/sec  and  the  basement 
velocity  is  4  km/sec.  As  in  the  Yucca  Flats  study,  the  calculations 
include  two  internal  multiple  reflections  within  the  structure.  The 
absolute  levels  of  motion  and  variability  of  ground  shaking  caused  by 
the  basin  structure  are  evidenced  in  the  displayed  synthetic 
seismograms  superimposed  at  seven  receiver  locations  across  the  basin. 
Note  the  long  duration  and  large  amplitudes  of  the  simulated  records 
along  the  left  half  of  the  basin  in  the  direction  of  the  incoming  energy. 
Amplifications  by  the  basin  structure  caused  by  focusing  is  responsible 
for  the  amplitude  variations  of  about  a  factor  of  three. 
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Figure  7.  Waveforms  and  relative  amplitudes  for  a  wave  incident  at  a 
shallow  angle  from  the  right  (arrow).  The  synthetic  time 
histories  have  been  calculated  with  the  Kirchhoff  technique 
and  includes  two  multiples  within  the  basin.  Note  the 
factor  of  3  increase  in  amplitude  caused  by  structure. 
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4.„  SOLUTION  USING  RIGOROUS  ITERATIVE  BIE  APPROACH 
In  this  chapter,  a  rigorous  iterative  approach  is  used  to  solve  the 
general  BIE  problem  formulated  in  Chapter  2.  The  solution  technique  is 
described  in  Section  4.1.  Included  is  a  discussion  of  how  the 
singularities  are  handled  in  the  general  formulation  which  leads  to  an 
iterative  treatment  based  on  perturbing  the  model  from  a  flat  layer 
point  of  view.  The  matrix  equation  becomes  well  conditioned  and 
numerically  tractable  after  the  singularities  have  been  removed  and  all 
matrix  inversion  operations  are  reduced  to  simpler  deconvolutional 
operations  by  virtue  of  the  iterative  solution  approach.  The  solution 
provides  kinematic  and  dynamic  effects  equally  accurately  at  low  and 
high  frequencies  and  handles  incident  waves  beyond  critical  angle  and 
even  layer  pinchouts. 

Sample  applications  are  presented  in  Section  4.2  for  the  simple  cases  of 
waves  impinging  on  an  irregular  free-surface  of  a  semi-infinite  half¬ 
space  and  transmission  through  an  irregular  interface  for  a  single  layer 
overlying  a  semi-infinite  half-space.  The  convergence  of  the  BIE 
solution  is  demonstrated  and  comparisons  are  shown  with  the 
approximate  Kirchhoff  solution.  The  deficiencies  of  the  Kirchhoff 
approximation  are  easily  identified  from  these  simple  models  and  the 
zero  order  iteration  is  shown  to  provide  results  more  accurate  than  the 
approximate  Kirchhoff  solution.  Applications  to  more  complicated  models 
are  discussed  in  Chapter  5. 

4.1  Methodology 

In  the  iterative  frequency  domain  treatment  of  the  general  BIE  problem 
presented  in  Chapter  2,  the  interaction  of  the  propagating  wavefield 
with  the  boundaries  is  satisfied  globally  including  rigorous  interaction 
between  neighboring  points  on  a  single  boundary.  The  first  step  in 
deriving  the  iterative  BIE  solution  is  to  remove  the  singularities  from 
the  Green's  functions  in  Eq.  (5).  These  singularities  occur  when  a 
quadrature  point  along  boundary  £  coincides  with  a  sample  point  along 
boundary  £  and  are  evidenced  along  the  diagonal  elements  of  sub¬ 
matrices  [G1^]  and  [H^  for  £=2, 3, . . . ,  N+1  and  n=£-1,£  (eg.,  when 


SGI -R-83-083 


39 


sample  point  x. 


coincides  with  quadrature  point  in 


the  second 


integral  of  Eq.  (4a)  and  the  first  integral  of  Eqs.  (2)  and  (4b)).  This 
is  typically  the  case  in  the  general  BIE  formulation  because  it  is 
desirable  to  be  able  to  use  the  same  sample  points  as  quadrature  points 
in  propagating  the  boundary  values  from  interface  £  to  interfaces  £-1 
and  £+1.  Singularities  can  also  occur  when  two  interfaces  pinchout  so 
that  a  quadrature  point  on  an  adjacent  boundary  (£+1  or  £-1)  coincides 
with  a  sample  point  along  boundary  £.  If  pinchouts  exist  in  the  model, 
these  singularities  would  be  evidenced  in  certain  elements  of 
submatrices  [G^+1  and  [H^+1  ^  or  [G^_n  £]  and  [H^_1  £]  (eg., 
when  sample  point  x^  coincides  with  quadrature  points  y  ^  or  y^.-j  in 
the  first  integral  of  Eq.  (4a)  or  in  the  second  integral  of  Eqs.  (2)  or 
(4b),  depending  on  which  two  layers  are  pinching  out).  Even  if  the 
sample  points  were  slightly  offset  from  the  quadrature  points,  the 
matrix  would  still  be  too  ill-conditioned  to  be  numerically  tractable. 


Basically  there  are  three  techniques  for  dealing  with  these  singularities 
with  various  degrees  of  generality  and  effectiveness.  One  method 
would  be  to  deal  with  the  singularities  directly  by  invoking  various 
smoothness  criteria  and  taking  principal  values.  This  has  the  distinct 
disadvantages  of  reducing  flexibility  in  the  formulation  and  of  being 
unstable  and  inefficient  during  applications.  A  second  method  would  be 
to  express  the  unknown  displacement  and  traction  boundary  values  in 
terms  of  integral  representations  of  unknown  forces  located  on  deformed 
surfaces  interior  to  the  actual  bounding  surfaces  and  then  to  solve  for 
the  set  of  unknown  forces  that  simultaneously  satisfies  all  the  boundary 
and  continuity  conditions  (see,  for  example,  Ohsaki,  1973  and  Part  II  of 
Apsel,  1979).  In  general,  this  treatment  is  quite  effective  because  the 
singularities  are  avoided  altogether  except  for  layer  pinchouts  (or  very 
thin  layers)  because  some  quadrature  points  on  the  deformed  surface 
would  still  coincide  some  of  the  sample  points  (or  be  so  close  that  the 
system  of  equations  becomes  ill-conditioned).  If  it  weren't  for  the 
possibility  of  layer  pinchouts,  the  equivalent  force  treatment  would  be 
useful  because  of  the  simpler  formulation  and  the  elimination  of  the  need 
to  evaluate  and  deal  with  the  spatial  derivatives  of  the  stress  Green's 
function  components  in  the  traction  integral  representations. 
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The  third  method,  which  is  the  preferred  choice,  involves  analytically 
subtracting  off  the  singularities  from  the  Green's  functions  using  the 
flat  layer  Green's  functions  and  iteratively  solving  the  complete  BIE 
problem  with  the  flat  layer  solution  as  a  starting  solution  on  the  first 
iteration.  This  has  the  distinct  advantages  of  removing  the 
singularities  under  ail  situations,  properly  conditioning  the  dense  BIE 
matrix  and  replacing  the  matrix  inversion  operations  involving  the  dense 
BIE  matrix  with  simple  deconvolution  operations  that  are  handled  quite 
efficiently  using  Fourier  transform  techniques. 

Another  important  advantage  of  using  this  approach  to  remove  the 
singularities  is  the  simple  provision  for  dealing  with  potentially  spurious 
edge  effects  arising  from  the  necessarily  finite  horizontal  extremes  of 
the  model  when  carrying  out  the  successive  iterations  on  the  solution. 
For  layer  boundaries  converging  to  the  depth  of  the  artificial  flat  layer, 
the  integrals  along  those  extremities  go  identically  to  zero  by  virtue  of 
the  perturbation  formulation.  For  those  boundaries  converging  to  a 
constant  depth  different  than  that  of  the  flat  layer,  the  first  order 
edge  effects  are  still  removed  leaving  higher  order  effects  governed  by 
the  ratio  of  the  horizontal  distance  from  sample  point  to  quadrature 
point  divided  by  the  vertical  offset  from  the  flat  layer  (indicating  that 
the  higher  order  edge  effects  can  be  removed  either  by  integrating 
along  this  extremity  until  the  offset  ratio  is  sufficiently  small  or  by 
tapering  the  diminishing  Green's  function  perturbations).  For  those 
boundaries  deviating  from  the  artificial  flat  layer  at  the  extremities  of 
the  model,  tapering  will  be  necessary  on  the  nonvanishing  Green's 
function  perturbations  unless  the  model  can  be  extended  until  the  layer 
boundaries  eventually  approach  zero  dip.  For  the  zero  order  iteration, 
tapering  will  typically  be  required  at  the  extremities  of  the  artificial  flat 
boundaries  to  obtain  the  flat  layer  solution  (unless  geometrical 
spreading  or  material  attenuation  behavior  are  sufficient  to  truncate  the 
integrals  without  tapering).  In  any  case,  it  is  quite  straightforward  to 
avoid  spurious  edge  effects  in  Sierra  Geophysics'  iterative  BIE  solution. 
This  is  in  direct  contrast  to  finite  difference  of  finite  element 
techniques  or  other  BIE  solution  techniques  where  either  appropriate 
nonreflecting  boundary  conditions  must  be  used  (such  as  in  Scheuster's 
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(1983)  BIE  solution)  or  else  the  model  grid  must  be  extended  far 

beyond  the  region  of  interest  (such  as  in  Cole's  (1980)  time  domain  BIE 
solution)  so  that  the  spurious  edge  reflections  arrive  sufficiently  late 
not  to  contaminate  the  signal  of  interest. 

The  iterative  solution  technique  at  a  given  frequency  for  the  general 
BIE  formulation  in  Eq.  (5)  is  derived  by  first  subtracting  off  the 
singularities.  This  is  effectively  accomplished  by  reformulating  the 
general  BIE  problem  in  terms  of  a  perturbation  problem  with  respect  to 
flat  layers.  The  matrix  equation  for  the  flat  layer  problem  in  which  the 
flat  layers  are  passed  through  the  reference  planes  of  the  original 
irregular  layers  is  written  as 

UHU0}  =  (H0HUo}  -  [Go]{To}  +  {FJ  (12) 

in  which  {UQ},  {Tq},  ^Fo^'  ]  and  [HQ]  have  the  analogous 

definitions  as  in  Eq.  (5)  for  the  flat  layer  problem  instead  of  the 

irregular  layered  problem.  In  the  limit  as  sample  point  and  quadrature 
point  coincide,  the  singular  elements  of  [Gq]  and  [Hq]  exactly  approach 
the  singular  elements  of  [G]  and  [H],  respectively.  Therefore,  if  the 
original  BIE  problem  could  be  recast  in  terms  of  the  perturbation 
matrices  [G-Gq]  and  [H-Hq],  formed  by  analytically  subtracting  the 
nontrivial  elements  of  [Gq]  and  [ Hq]  from  the  corresponding  elements  of 
[G]  and  [H],  then  all  singularities  would  be  automatically  handled. 

As  an  example,  consider  the  form  of  the  singularities  arising  from 
interaction  of  the  wavefield  with  a  boundary  in  2-D  acoustic  structures. 
The  horizontal  sample  position  on  the  irregular  boundary  is  defined  by 
the  coordinates  (x,z),  with  z  being  a  function  of  x.  The  artificial  flat 

layer  is  defined  at  a  constant  depth  of  z  =  zq  and  integration  is  carried 

out  along  the  boundary  over  all  quadrature  points  xQ.  The  singularity 
arises  in  the  imaginary  part  of  the  Hankel  functions  of  order  zero  and 
one  when  quadrature  point  xq  approaches  sample  point  x  in  [G]. 
Then,  when  the  perturbation  matrix  [G-Gq]  and  [H-HqJ  are  formed,  the 
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difference  of  the  imaginary  parts  of  the  singular  elements  (irregular 
operator  minus  flat  operator)  has  the  following  form  including  the 
quadrature  coefficients: 

{  Ym[k7(x-xo)2+(z-zo)2]7l+(dz/dx)2  -  Ym[ k |x-xq1  ]  }  dxQ  ,  (m=0,l)  (13) 

in  which  Y  (r)  is  the  Bessel  function  of  the  second  kind  of  order  m 
and  k  is  the  horizontal  wavenumber  equal  to  frequency  divided  by 
acoustic  velocity.  Using  the  fact  that  Yg(r)  ~  (2/n)  ln(r)  and  that 
Y^(r)  ~  -(2/n)/r  in  the  limit  as  r  goes  to  zero,  it  is  straightforward  to 
prove  (using  I'Hospital's  Rule)  that  the  difference  of  the  two  terms  in 
Eq.  (13)  is  identically  zero  in  the  limit  as  x  approaches  xq.  The  entire 
row  of  terms  appearing  in  matrices  [G-Gq]  and  [H-H  ]  for  other  values 
of  xq  in  Eq.  (13)  are  in  fact  better  conditioned.  The  behavior  at  the 
extremities  of  the  boundary  is  also  evident  from  Eq.  (13).  For 

irregular  boundaries  converging  to  the  artificial  flat  layer,  coordinate  z 
approaches  the  constant  zq  and  the  difference  in  Eq.  (13)  is  identically 
zero;  for  boundaries  converging  to  a  different  constant  depth,  Zy  the 
difference  approaches  zero  as  the  absolute  value  of  the  ratio 
(x-xo)/(z^-zQ)  increases. 

The  iterative  BIE  solution  technique  proceeds  by  rewriting  Eq.  (5) 
without  change  by  adding  and  subtracting  the  flat  layer  terms: 


[I]{U+Uo}  =  [H+Hq] {U+Uq}  -  (G+GoJ{?+To}  +  {F+Fq}  (14) 

in  which 

{U}  =  {U-Uo>  ,  {T}  =  {T-Tq}  (15) 

are  the  unknown  perturbed  displacement  and  traction  boundary  value 
vectors,  respectively; 

[G]  =  IG-GJ,  [H]  =  [H-Ho]  (16) 
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are  the  perturbed  displacement  and  traction  Green's  function  matrices; 
and 


{F}  =  {F-Fq} 


(17) 


is  the  perturbed  direct  source  field  vector.  Eq.  (12)  is  then  used 
directly  to  simplify  Eq.  (14),  which  after  rearranging,  takes  the  form 


[I] {U}  -  [Ho]{U}  +  [GoHT}  =  lH]{U+Uo}  -  [G]{f+To}  +  {F},  (18) 

in  which  all  matrices  are  nonsingular  and  well  behaved:  matrix  [I] 
consists  of  identity  and  null  submatrices;  the  rows  of  matrices  [Gq]  and 
[Hq]  are  convolutional  operators  on  vectors  {T}  and  {U}  for  constant 
material  properties  in  each  layer;  and  the  rows  of  matrices  [H]  and  [ G ] 
are  small  and  smoothly  varying  for  structures  not  deviating  too  widely 
from  the  reference  flat  layers. 

Next,  Eq.  (18)  is  rewritten  using  more  compact  ontation: 


[C] {X}  =  [Al{X+Xo}  +  {F}  (19) 

in  which  vector  {X}  contains  all  the  unknown  perturbed  boundary 
values  along  all  layer  boundaries: 


x  =  {u1},{u2},{t2},{u3},{t3} . {Un+iU?N+1}  (20) 

where  subvectors  {U£}  contains  the  3q£  unknown  perturbed  dis¬ 
placement  boundary  values  (u-uQ)  along  interface  £  (£=1 ,2, . . .  ,N+1 )  and 
subvectors  contain  the  3q^  unknown  perturbed  traction  boundary 

values  (T-T  )  along  interface  £  (£=2,3, . . .  ,N+1 ).  Vector  {XQ}  is  de¬ 
fined  analogously  to  Eq.  (20)  for  the  flat  layer  displacement  and  trac¬ 
tion  boundary  values,  so  that  the  actual  boundary  values  are  given  by 
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{X}  =  {X}  +  {Xo}.  (21) 

The  square  matric  [A]  contains  the  relatively  small  and  well  behaved 
submatrices  from  [H]  and  [G].  The  square  matrix  [C]  contains  the 
identity  submatrices  from  [I]  and  the  convolutional  submatrices  from 
IG  ]  and  [Hq]  .  Therewith  Eq.  (19)  can  be  written  in  its  final  form, 
which  is  ideally  suited  for  an  iterative  solution  approach: 

{X}n  =  {X}0  +  [p]{Xn_i  +  Xo}  >  n  =  1,2,...  (22) 

in  which  the  zero  order  iterative  solution  {X}g  is  given  by 

{X}Q  =  l C ] ~ 1 {F }  ;  (23a) 

the  flat  layer  solution  {XQ}  is  given  by 

{Xq}  =  [C]“1{Fo>  ;  (23b) 

and  the  recursive  matrix  operator  [P]  is  given  by 

[P]  =  I C ] ~ 1 [A] .  (23c) 

Note  that  {X}q,  {Xq}  and  [P 1  are  all  independent  of  iteration  number  n 
and  are  efficiently  calculated  all  at  once  in  the  Fourier  transform  domain 
as  the  inverse  operation  involving  matrix  [C]  decouples  into  a  simpler 
deconvolutional  problem. 

The  nth  iteration  estimates  to  the  perturbed  boundary  values  are 
obtained  by  recycling  the  values  from  iteration  number  n-1  through  the 
right-hand  side  of  Eq.  (22)  to  obtain  an  improved  solution  to  the  un- 
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known  perturbed  boundary  values.  The  recursive  process  in  Eq.  (22) 
is  numerically  stable  and  convergent.  The  convergence  proof  would 
require  a  detailed  discussion  of  the  norm  of  matrix  [P]  in  Eq.  (23c) 
and  will  not  be  presented  here.  However,  in  all  cases  where  the 
technique  has  been  applied,  the  successive  iterations  have  exhibited 
convergent  behavior.  As  the  number  n  of  successive  iterations  is 
increased,  the  amount  of  high  order  interaction  of  the  wavefield  with 
the  boundaries  is  also  increased.  Also,  the  zero  order  iteration  is 
often  considerably  more  accurate  than  the  approximate  Kirchhoff 
solution  of  the  previous  chapter  and  is  useful  in  itself. 

Once  the  desired  number  of  iterations  has  been  carried  out,  the  final 
solution  for  the  unperturbed  boundary  values  is  obtained  from  Eq.  (21) 
in  which  the  final  solution  to  the  perturbed  boundary  values  {X}  is 
given  by  the  nth  order  iteration  in  Eq.  (22)  and  the  flat  layer  solution 
{XQ}  is  given  by  Eq.  (23b).  Then,  analogous  to  the  final  propagation 
step  in  the  Kirchhoff  solution  technique,  the  displacement  field  at  any 
additional  locations  within  any  layer  is  obtained  at  this  frequency  by 
introducing  the  boundary  values  from  Eq.  (21)  into  the  representation 
theorem  of  Eq.  (1).  Time  domain  results,  if  desired,  would  be  obtained 
through  Fourier  synthesis  by  calculating  the  iterative  BIE  solution  at  a 
discrete  set  of  frequency  values. 


4.2  Results 

An  iterative  frequency  domain  BIE  algorithm  has  been  developed  for 
treating  three  simple  cases  using  the  methodology  outlined  in 
Section  4.1.  The  three  cases  include  the  following  two-dimensional 
acoustic  problems:  (1)  wave  propagation  in  a  semi-infinite  space  with 
an  irregular  free  surface;  (2)  wave  propagation  through  an  irregular 
interface  in  an  infinite  space;  and  (3)  wave  propagation  in  a  layer 
overlying  a  semi-infinite  space  with  an  irregular  interface  and  a  flat 
free  surface.  The  results  for  these  simple  cases  illustrate  the 
effectiveness  of  the  iterative  BIE  technique  and  are  presented  in  this 
section.  Extensions  to  a  more  general  iterative  BIE  algorithm  are 
discussed  in  Chapter  5. 
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A  comparison  between  the  Kirchhoff  solution  (from  Section  3.1)  and  the 
iterative  BIE  solution  (from  Section  4.1)  is  presented  in  Figures  8  and 
9  for  the  case  of  a  plane  monochromatic  wave  in  a  semi-infinite  space 
normally  incident  to  an  irregular  free  surface  with  a  simple  bump.  The 
power  density  of  the  scattered  field  is  shown  in  Figure  8  at  a 
frequency  corresponding  to  a  horizontal  wavenumber  (to/a)  of  1.5  for 
the  half-space.  The  Kirchhoff  solution  is  shown  in  the  upper  left  plot 
with  the  zero  order,  first  order  and  successive  iterations  of  the  BIE 
solution  shown  in  the  center  left,  lower  left  and  right-column  plots, 
respectively.  In  each  of  the  six  plots,  the  response  is  shown  as  a 
function  of  increasing  horizontal  position  and  increasing  vertical  depth 
into  the  half-space  with  the  amplitudes  color  scaled  between  dark  red 
for  the  maximum  value  through  yellow,  green  and  light  blue  for  the 
intermediate  values,  down  to  dark  blue  for  the  minimum  value  (zero  for 
the  power  density  plots).  It  is  seen  that  the  Kirchhoff  solution  fails  to 
reproduce  the  diffractive  effects  arising  along  the  incline  of  the  surface 
bump.  By  the  second  iteration  of  the  iterative  BIE  soluton,  these 
diffractive  effects  are  nicely  produced  and  convergence  is  verified  by 
comparison  to  the  third  and  fourth  iteration  plots.  Even  the  zero  order 
BIE  solution  is  considerably  more  accurate  than  the  Kirchhoff  solution. 

The  corresponding  time  domain  comparisons  are  shown  in  Figure  9.  For 
the  time  domain  plots,  the  vertical  axis  represents  the  horizontal 
receiver  positions  as  a  function  of  increasing  distance  along  a  flat  line 
several  grid  dimensions  below  the  reference  plane  of  the  irregular  free 
surface.  The  color  coding  convention  is  the  same  as  in  "igure  8  with 
red  being  the  most  positive  (at  the  focus  of  scattered  energy)  and  dark 
blue  being  the  most  negative  (along  the  reflected  power  plane  wave). 
The  response  has  been  convolved  with  a  Gaussian-shaped  wavelet.  The 
largest  deficiencies  evidenced  in  the  Kirchhoff  solution  in  the  upper  left 
plot  are  the  inaccurate  amplitudes  and  inaccurate  travel-time  delays  for 
the  later  diffractive  arrivals,  substantiating  the  indicated  deficiencies  at 
the  particular  frequency  used  in  Figure  8.  The  excellent  convergence 
of  the  BIE  solution  in  the  time  domain  verifies  that  the  technique  works 
equally  well  at  low  and  high  frequency.  Also,  there  are  no  edge 
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Figure  8.  Comparison  of  the  power  density  of  the  scattered  field  due 
to  a  vertically  incident  plane  wave  impinging  on  the  free 
surface  with  a  bump  between  the  Kirchhoff  solution  and  the 
various  iterations  of  the  more  rigorous  BIE  solution.  All 
plots  are  scaled  to  absolute  range  (0.1,  3.0). 
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Figure  9.  Corresponding  comparison  to  Figure  8  in  the  time  domain 
between  the  Kirchhoff  solution  and  the  iterative  BIE 
solutions  along  a  line  of  receivers  located  several  grid 
dimensions  below  the  free  surface.  All  plots  are  scaled  to 
absolute  range  (-1.0,  0.5). 
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reflections  contaminating  the  signal  verifying  that  the  technique 

effectively  avoids  these  spurious  arrivals.  Again,  even  the  zero  order 
BIE  solution  is  considerably  more  accurate  than  the  Kirchhoff  solution. 

The  next  three  figures  show  the  sensitivity  of  the  iterated  BIE  solution 
to  various  model  parameters  for  the  case  of  a  simple  bump  on  the  free 
surface  of  a  half-space.  In  Figure  10,  the  power  density  of  the 

scattered  field  due  to  a  normally  incident  plane  wave  is  compared  at  six 
frequencies  from  low  frequency  (upper  left  plot  with  k=0.1)  to  high 
frequency  (lower  right  plot  with  k=2.1)  with  each  plot  self-scaled.  The 
most  obvious  trend  to  notice  in  these  six  plots  is  the  increasingly 
focussed  reflections  and  diffractions  from  the  surface  bumps  as  the 
frequency  is  increased.  At  low  frequency,  the  wavelength  is 

approaching  or  exceeding  the  dimensions  of  the  surface  bump,  reducing 

the  diffractive  effects  and  causing  the  bump  to  act  like  a  point  source 
emitter.  In  Figures  11  and  12,  the  iterated  BIE  solution  due  to  a 
normally  incident  plane  wave  is  compared  for  different  size  surface 
bumps  in  the  frequency  domain  and  time  domain,  respectively.  The 
power  density  of  the  scattered  field  at  an  intermediate  frequency 
corresponding  to  k=1.5  is  shown  for  no  bump  in  the  upper  left  hand 
plot  to  the  largest  bump  in  the  lower  right  hand  plot,  with  all  plots 
displayed  using  the  same  absolute  color  scaling.  As  expected,  there  is 
no  scattered  field  for  the  flat  free  surface  with  increasingly  large-' 
scattered  energy  and  larger  dispersion  of  diffractive  energy  as  the 
bump  is  increased  in  size.  The  analogous  results  in  the  time  domain 
including  the  reflected  wave  are  shown  in  Figure  12  for  a  line  of 
receivers  located  several  grid  dimensions  below  the  free  surface.  All 
six  plots  are  displayed  using  the  same  absolute  color  scaling  with  the 
most  negative  arrivals  (dark  blue)  corresponding  to  the  reflected  plane 
wave  and  the  most  positive  arrivals  (red)  corresponding  to  the  maximum 
scattered  field. 

The  next  four  figures  show  the  sensitivity  of  the  iterated  BIE  solution 
to  various  model  parameters  for  a  more  arbitrary  irregular  free  surface. 
In  Figures  13  and  14,  the  scattered  power  density  field  is  shown  for 
vertically  incident  and  obliquely  incident  (30  degrees  from  vertical) 
plane  waves,  respectively.  In  both  figures,  the  field  is  compared  at 
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Figure  10.  Power  density  of  the  scattered  field  due  to  a  vertically 
incident  plane  wave  impinging  on  a  free  surface  v\ith  a 
bump  calculated  with  the  iterative  BIE  technique  at  six 
different  frequencies.  All  plots  are  self-scaled. 
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Figure  11.  Comparison  of  the  power  density  of  the  scattered  field  for 
the  free  surface  with  a  bump  as  a  function  of  the  amplitude 
of  the  bump.  All  plots  are  scaled  to  absolute  range 
(0.1,  3.0). 
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Figure  12.  Corresponding  results  to  Figure  11  in  the  time  domain  for  a 
line  of  receivers  located  several  grid  dimensions  below  the 
free  surface.  All  plots  are  scaled  to  absolute  range  (-1.0, 
0.5). 
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Figure  13.  Power  density  of  the  scattered  field  due  to  a  vertically 
incident  plane  wave  impinging  on  a  highly  irregular  free 
surface  calculated  with  the  iterative  BIE  technique  at  four 
different  frequencies.  All  plots  are  self-scaled. 
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Figure  14.  Analogous  results  to  Figure  13  for  a  plane  wave  incident  at 
30  degrees  from  vertical.  All  plots  are  self-scaled. 
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four  different  frequencies  corresponding  to  wavenumbers  of  0.1,  0.5, 
1.0  and  2.0  in  the  upper  left,  lower  left,  upper  right  and  lower  right 
plots,  respectively.  As  anticipated  from  the  analogous  sensitivity  study 
in  Figure  10  for  the  case  of  a  simple  bump  irregularity,  increasingly 
focussed  reflections  and  diffractions  are  observed  from  the  surface 

irregularities  as  the  frequency  is  increased.  Again,  at  low  frequency, 
the  wavelengths  are  approaching  or  exceeding  the  dimensions  of  the 
surface  irregularities,  which  reduces  the  diffractive  effects  and  causes 
the  irregularities  to  resemble  local  point  source  emitters.  Time  sections 
for  both  angles  of  incidence  considered  in  Figures  13  and  14  are 

displayed  in  Figure  15  for  a  line  of  receivers  located  just  below  the  free 
surface.  The  normal  incidence  section  is  on  the  left  and  the  oblique 

incidence  section  is  on  the  right  and  are  in  the  same  format  as  the 

plots  of  Figures  9  and  12.  As  anticipated  from  the  power  density  fields 
at  the  sampling  of  frequencies  shown  in  Figures  13  and  14,  all 
diffractive  effects  are  accurately  simulated  in  the  iterative  BIE  solution 
and  there  are  no  spurious  edge  effects. 

In  Figure  16,  the  sensitivity  of  the  time  section  for  the  obliquely 
incident  plane  wave  case  of  Figure  15  is  examined  with  respect  to 
amplitude  of  the  surficial  irregularities.  The  four  time  sections 
represent  the  iterated  BIE  response  compared  with  surficial 
irregularities  of  four  different  amplitudes,  respectively.  All  four  plots 
are  displayed  using  the  same  absolute  color  scaling.  In  the  lower  left 
plot,  the  surficial  irregularities  are  sufficiently  small  relative  to  the 
shortest  wavelengths  considered  in  the  calculation  that  the  scattered 
energy  is  negligible  compared  to  the  primary  reflected  wave.  As  the 
surficial  irregularites  are  increased,  the  scattered  energy  increases  and 
the  diffracted  waves  are  seen  to  be  as  large  in  absolute  value  as  the 
primary  reflected  wave. 

The  case  of  propagation  through  an  interface  is  considered  in  Figures 
17  through  20.  The  first  two  of  these  figures  demonstrate  the  accuracy 
of  the  iterative  BIE  technique  for  this  case.  In  Figure  17,  the  power 
density  of  the  scattered  transmitted  field  through  the  interface  is 
compared  using  the  zero  order  BIE  solution  for  an  impedance  contrast 
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Figure  15.  Iterative  BIE  time  domain  solution  corresponding  to  Figures 
13  and  14  for  a  line  of  receivers  located  several  grid 
dimensions  below  the  free  surface.  Both  plots  are 
self-scaled . 
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Figure  16.  Comparison  of  iterated  time  domain  solution  for  the  highly 
irregular  free  surface  as  a  function  of  the  amplitude  of  the 
irregularities .  All  plots  are  scaled  to  absolute  range  (-1.0, 
0.5). 
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Figure  17.  Comparison  of  power  density  of  the  scattered  field  due  to 
propagation  through  an  interface  with  and  without  an 
impedance  contrast  using  the  zero  order  BIE  solution.  Both 
plots  are  scaled  to  absolute  range  (0.0,  3.5). 
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Figure  18.  Magnification  of  power  density  of  the  scattered  field  for  the 
zero  impedance  contrast  case  from  Figure  17  as  a  function 
of  iteration  number,  showing  convergence  to  null  solution. 
Both  plots  are  scaled  to  absolute  range  (0.0,  0.2). 
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Figure  19.  Power  density  of  the  scattered  field  due  to  propagation 
through  an  interface  as  a  function  of  angle  of  plane  wave 
emergence  (critical  angle  at  30  degrees).  All  plots  are 
self-scaled. 
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Figure  20.  Power  density  of  the  scattered  field  due  to  propagation 
through  an  interface  at  intermediate  and  high  frequencies. 
All  plots  are  self-scaled. 
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of  1.5  in  the  left  plot  and  no  impedance  contrast  in  the  right  plot. 
Both  plots  are  displayed  on  the  same  absolute  color  scale  and  as 
expected,  no  transmitted  scattered  energy  is  displayed  for  zero 
impedance  contrast.  The  field  for  the  zero  impedance  contrast  case  is 
magnified  in  the  left  plot  of  Figure  18  showing  the  numerically  null 
solution.  While  the  field  is  negligible  compared  to  the  finite  impedance 
contrast  case  in  Figure  17,  the  numerical  solution  is  never  exact. 
Fortunately,  however,  the  iterated  BIE  solution  converges  to  the  exact 
solution  as  more  and  more  iterations  are  performed  as  demonstrated  in 
the  right  plot  of  Figure  18  with  the  next  order  iterated  solution 
converging  an  order  of  magnitude  closer  to  the  exact  solution. 

In  Figures  19  and  20,  the  acoustic  velocity  above  the  interface  is  twice 
the  velocity  below  the  interface,  so  that  the  critical  angle  of  emergence 
is  30  degrees  from  vertical.  The  power  density  of  the  scattered 
transmitted  field  at  an  intermediate  frequency  is  shown  in  Figure  19  as 
a  function  of  emergence  angle.  As  demonstrated  by  the  sequence  of 
plots  for  a  planar  interface  from  pre-critical  emergence  in  the  left 
column  of  plots  to  post-critical  emergence  in  the  upper  right  plot,  the 
iterative  BIE  approach  provides  the  correct  results  under  all  cases.  In 
particular,  the  critically  refracted  waves  in  the  upper  right  plot 
properly  decay  as  a  function  of  distance  from  the  boundary  (note  that 
these  waves  can  not  be  simulated  using  the  approximate  Kirchhoff 
solution  technique).  The  post-critical  case  is  repeated  in  the  lower 
right  plot  for  a  simple  non-planar  interface.  The  important  feature  to 
notice  is  how  the  bump  acts  like  a  point  source  and  radiates  a 
significant  amount  of  energy  in  the  general  direction  of  the  incoming 
wave  that  would  otherwise  be  absent  as  shown  in  the  upper  right  plot. 
In  the  case  of  a  multilayered  half-space,  irregularities  along  interfaces 
will  give  rise  to  late  arriving  head  wave  energy  at  the  free  surface  that 
would  otherwise  not  be  evidenced  with  planar  interfaces.  In  Figure  20, 
the  power  density  of  the  scattered  transmitted  field  due  to  normally 
incident  propagation  through  simple  non-planar  interfaces  is  compared 
at  intermediate  and  high  frequency  (three  times  intermediate  value). 
The  comparison  is  performed  for  irregular  interfaces  consisting  of  a 
single  bump  and  a  double  bump  in  the  left  and  right  columns, 
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respectively.  The  important  feature  to  notice  here  is  the  significant 
beaming  of  the  scattered  field  at  high  frequency  (lower  row  of  plots)  in 
contrast  to  the  simple  point-source  like  emissions  from  the  irregularity 
at  low  frequency. 

The  case  of  propagation  through  an  irregular  interface  with  a  flat  free 
surface  overhead  is  considered  in  Figure  21.  The  acoustic  velocity  in 
the  layer  is  half  the  value  of  the  underlying  half-space.  The  power 
density  of  the  scattered  field  due  to  a  vertically  incident  plane  wave  is 
compared  at  intermediate  and  high  frequencies  (twice  the  intermediate 
value)  in  the  left  and  right  plots,  respectively,  for  a  relatively  thick 
layer.  As  in  Figure  20,  much  more  focussing  is  experienced  at  high 
frequency  (right  plot)  and  surface  waves  can  be  observed  along  the 
free  surface  at  lower  frequency  (left  plot)  which  decay  exponentially  as 
a  function  of  depth  for  the  thick  layer.  The  center  plot  is  the  same  as 
the  left  plot  except  that  the  layer  thickness  is  considerably  reduced 
giving  rise  to  a  large  trapped  surface  wave  propagating  bilaterally  away 
from  the  irregularity  for  the  normally  incident  plane  wave. 
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Figure  21.  Power  density  of  the  scattered  field  due  to 
propagation  through  an  irregular  interface  with  a 
flat  free  surface  overhead.  All  plots  are 
self-scaled. 


SGI -R -83-083 


65 


5.0  RECOMMENDED  EXTENSIONS  OF  THE  ITERATIVE  BIE  ALGORITHM 
An  iterative  BIE  algorithm  has  been  developed  for  acoustic,  two- 
dimensional  structures  having  no  more  than  a  single  irregular  layer 
overlying  a  semi-infinite  half  space.  The  simplicity  of  the  algorithm  has 
allowed  efficient  testing  of  the  methodology  and  refinement  of  the 
iterative  procedure.  Results  for  the  simple  cases  supported  under  this 
algorithm  were  presented  in  Section  4.2  and  are  quite  useful  in 
providing  new  insights  into  wave  propagation  problems  involving 
irregular  interfaces.  Now  that  the  method  has  been  shown  to  be 
reliable,  complete  and  highly  efficient  (relative  to  finite  difference  type 
techniques),  a  more  general  iterative  BIE  algorithm  is  needed  to  study 
more  complex  wave  propagation  problems. 

The  formulation  in  Section  4.1  is  completely  general  for  an  irregular, 
multilayered,  elastic,  three-dimensional  half-space  and  the  solution 
technique  is  outlined  in  detail.  The  most  straightforward  extension  of 
the  algorithm  is  to  the  case  of  a  multilayered,  two-dimensional  half¬ 
space  because  the  Green's  functions  are  the  same  as  for  the  single 
layer,  2-D  case  and  the  deconvolutional  form  of  the  BIE  matrix  equation 
can  still  be  handled  with  one-dimensional  Fourier  transform  techniques. 
However,  the  two-dimensional  class  of  problems  is  limiting  and  the  next 
recommended  extension  is  to  the  three-dimension  case.  The  extension 
to  the  full  multilayered,  3-D  case  is  straightforward  in  theory  with  the 
Green's  functions  simplifying  to  exponential  functions  from  Hankel 
functions  (refer  to  Appendix  A)  and  with  the  deconvolutional  form  of 
the  BIE  matrix  equations  handled  with  two-dimensional  Fourier  transform 
techniques.  The  efficiency  of  the  3-D  code  will  depend  very  strongly 
on  optimized  handling  of  the  large  core  requirements  and  spatial 
gridding  requirements  as  a  function  of  frequency.  It  is  anticipated 
that  with  sufficient  effort  on  developmental  efficiency,  the  code  will 
execute  up  to  an  order  of  magnitude  faster  than  finite  difference  type 
codes. 

It  is  recommended  that  the  matrix  formalism  be  retained  in  all 
extensions  of  the  algorithm  to  permit  optimally  efficient  execution 
through  use  of  array  processing  hardware  if  desired.  Also,  the  elastic 
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case  should  be  explicitly  supported  as  the  form  of  algorithms  is  the 
same  with  just  more  bookkeeping  required  for  the  expanded  set  of 
nontrivial  components  at  each  node,  and  a  number  of  source  types 
should  be  explicitly  supported  including  point  forces,  explosions,  point 
dislocations,  and  finite  sources. 
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APPENDIX  A.  Infinite-space  Green's  functions 

Expressions  for  the  infinite-space  Green's  functions  are  given  in  this 
appendix  for  two  and  three  dimensional  wave  propagation  problems  in 
acoustic  and  elastic  material.  In  each  case,  the  displacement  and  stress 
fields  at  location  y  due  to  an  impulsive  point  source  at  location  x  are 
presented  as  dimensionless  quantities  per  unit  force  in  the  frequency 
domain  for  compatibility  with  the  frequency  domain  BIE  method. 

1.  Three-dimensional,  elastic  case 

The  dimensionless  displacement  field  per  unit  force  is  given  by 


(4nMR)G£i„(x,y;u>)  =  [x£m(l+3da)  - 


-ik„ 


'  l*£m(1+3V  ■  6*m(1+V]  e 


(A-l) 


in  which  G£m(x,y;w)  is  the  ^-component  of  the  complex  displacement 
field  at  location  y  and  frequency  u>  due  to  a  point  source  in  the  in¬ 
direction  at  location  x  with  £,m=l,2,3.  The  first  main  term  on  the 
right-hand  side  corresponds  to  the  compressional  wave  traveling  at  a 
velocity  a  and  the  second  main  term  corresponds  to  the  shear  wave 
traveling  at  a  velocity  p.  The  various  quantitives  appearing  in 
Eq.  A-1)  are  defined  by  the  following  relations: 


x  -  (Xj  ,x£  ,x^)  =  source  location;  y  =  (yj^.Y-j)  =  receiver  location; 


R  = 


=  I  y-x  I  =  ((yj-Xj)2  +  (y2-x2)2  + 

=  Vm  5  =  (VXn,)/R  ’  m  =  1>2>3  ; 

=  1  if  2=m;  6£m  =  0  if  ***'> 

a 2  =  (m+2X)/y  ;  P2  =  m/P  ;  Y  =  0/“ 

M,X  =  Lame  constants  ;  p  =  density  ; 

ka  =  u>/a  ;  kp  =  ui/p  ;  i  =  (-1)^  ; 

da  =  (ikaR)_1  { l+(ikaR)_1 }  ;  dp  =  (ikpR)_1  (l  +  dk^R)'1] 
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In  the  limit  as  source  and  receiver  point  coincide,  R  goes  to  zero 
causing  singularities  of  the  form  R  *,  j=1,2,3  in  G^.  Conversely,  in 
the  limit  as  R  goes  to  infinity,  the  near-field  terms  involving  d^  and  d^ 
vanish  leaving  a  geometrical  attenuation  rate  of  (1/R).  Material 
attenuation  may  be  introduced  by  designating  imaginary  parts  to  the 
velocities  a  and  (3  so  that  the  exponential  functions  provide  further 
decay  at  large  distances.  Also,  these  imaginary  parts  may  be  specified 
with  arbitrary  frequency  dependence. 


The  corresponding  dimensionless  orders  stress  field  per  unit  force  is 
given  by 


4nR  Wx’y;u)) 


+ 


< -  x0  [2y2(ik  R  +  6  +  15d  )] 
l  &nnl  '  a  a'J 

+  d£mnt(2Y  '1)ikaR  +  <4Y  "D  +  6V 

{  W2<V  +  6  *  15dp>l  -  <W2<‘*3V1 

-  +  3  ♦  6<y }  *  p 


(A-2) 


in  which  G^mn(x,y;ui)  is  the  £m-component  of  the  complex  stress  field 
at  location  y  and  frequency  u>  due  to  a  point  source  in  the  n-direction 
at  location  x  with  £,m,n  =  1,2,3.  Again,  the  first  and  second  main 
terms  on  the  right-hand  side  correspond  to  the  compressional  and  shear 
waves,  respectively.  Those  quantities  appearing  in  Eq.  (A-2)  that 
were  not  used  in  Eq.  (A-1)  are  defined  by  the  following  relations: 


£mn 


=  X.X  X 

£  m  n 


d .  _  6„  x 

£ran  =  £m  n 


The  singularities  in  the  stress  components  are  of  the  form  R 

j=1 ,2,3,4  as  R  goes  to  zero  which  is  one  order  higher  than  the 
corresponding  singularities  in  the  displacement  components.  However, 
the  far-field  geometrical  (and  material)  attenuation  behavior  is  the  same 
as  described  for  the  displacements. 
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The  £-component  of  the  complex  traction  field  at  location  y  is  then 
obtained  using  Cauchy's  formula  by  forming  the  dot  product  of  the 
stress  tensor  with  the  unit  normal  vector  v  at  that  point: 


H£n(x’y;U)) 


Gm£n(x-y;u0 


vJy) 


m=l ,2,3. 


(A-3) 


2.  Three-dimensional,  acoustic  case 

The  dimensionless  displacement  field  per  urit  force  is  given  directly 
from  Eq.  (A-1)  to  read: 


[47i(\+2M)R]G£m(x,y;u.)  =  [x£m(l+3da)  -  fi^dj  e  “  (A-4) 

with  the  corresponding  dimensionless  stress  field  per  unit  force  given 
directly  from  Eq.  (A -2)  to  read: 

+  6  4  I5V> 

*  dto»l<2^-1)ikaR  *  +  V 

4  lko,,!  <A-5) 

The  traction  field  may  then  be  evaluated  by  introducing  Eq.(A-5)  into 
Eq.  (A-4). 

3.  Two-dimensional  elastic  case 


The  Green's  functions  are  presented  for  the  two-dimensional  case  of 
antiplane  strain.  There  is  only  one  nonzero  displacement  component 
and  one  nonzero  traction  component  for  antiplane  strain.  Assuming  that 
*2  is  the  translational  direction,  the  the  dimensionless  displacement 
component  per  unit  force  is  given  by 
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4(jG22(X1  ,x;};y1  ,y3;u))  =  -iy2H^  ^(k^r)  -  iH^  ^(k^r) 


(A-6) 


in  which  H 


(2) 

'n 


is  the  Hankel  function  of  the  second  kind  of  order  n  and 

'z 


r  =  [(x^-y^)  +  (x3-y3)  ]  2.  The  corresponding  dimensionless  traction 

component  per  unit  force  is  given  by: 


(2) 

4rH22(x1,x3;y1,y3;uj)  =  (kar) 

(2\ 

+  ikp(y3-x3)H^  J(kpr) 


(A- 7 ) 


In  this  case,  the  singularities  arising  when  source  and  receiver  points 
coincide  (r=0),  occur  for  the  displacement  component  in  Eq.  (A-6)  in 
the  imaginary  part  of  the  zero  order  Hankel  function,  exhibiting 
logarithmic  behavior  (i.e.,  Bessel  function  of  the  second  kind 
Yq(z)  ~  (2/7i)ln(z)  as  z  approaches  zero).  For  the  stress  components 
in  Eq.  (A-7),  the  singularities  occur  in  the  imaginary  part  of  the 
Hankel  functions  of  order  1  (i.e.,  Bessel  function  Y^(z)  ~  -(2/n)/z  as  z 
approaches  zero).  The  cases  of  plane  strain  or  plane  stress  exhibit 
similar  singularity  behavior. 


4.  Two-dimensional  acoustic  case 

The  displacement  and  traction  components  are  given  directly  from 
Eqs.  (A-6)  and  (A-7)  to  read: 


4  ( \+2p)  G22  (x x ,  x2 ; y x ,  y3 ; ui) 
4rH22(x1,x3;y1,y3;ui) 


( 2\ 

'“o  (k»r) 

ik0(y3-*3)  HpOy) 


(A-8) 

(A-9) 
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